CMS 3D CMS Logo
Primary4DVertexValidation.cc
Go to the documentation of this file.
1 #include <numeric>
2 #include <string>
3 #include <vector>
4 #include <map>
5 #include <algorithm>
6 
7 #include "FWCore/Framework/interface/Frameworkfwd.h"
8 #include "FWCore/Framework/interface/ESHandle.h"
9 
10 #include "DataFormats/Common/interface/ValidHandle.h"
11 #include "FWCore/Framework/interface/Event.h"
12 #include "FWCore/Framework/interface/EventSetup.h"
13 #include "FWCore/Framework/interface/MakerMacros.h"
14 #include "FWCore/ParameterSet/interface/ParameterSet.h"
15 #include "FWCore/MessageLogger/interface/MessageLogger.h"
16 
17 #include "DataFormats/Math/interface/LorentzVector.h"
18 #include "DataFormats/Math/interface/Point3D.h"
19 
20 // reco track and vertex
21 #include "DataFormats/TrackReco/interface/Track.h"
22 #include "DataFormats/TrackReco/interface/TrackFwd.h"
23 #include "DataFormats/VertexReco/interface/Vertex.h"
24 #include "DataFormats/VertexReco/interface/VertexFwd.h"
25 #include "RecoVertex/VertexPrimitives/interface/TransientVertex.h"
26 
27 // TrackingParticle
28 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
29 #include "SimDataFormats/Associations/interface/TrackToTrackingParticleAssociator.h"
30 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h"
31 #include "SimDataFormats/TrackingAnalysis/interface/TrackingVertexContainer.h"
32 
33 // pile-up
34 #include "SimDataFormats/PileupSummaryInfo/interface/PileupSummaryInfo.h"
35 
36 // associator
37 #include "SimTracker/VertexAssociation/interface/calculateVertexSharedTracks.h"
38 
39 // vertexing
40 #include "RecoVertex/PrimaryVertexProducer/interface/TrackFilterForPVFinding.h"
41 
42 // simulated vertex
43 #include "SimDataFormats/Associations/interface/VertexToTrackingVertexAssociator.h"
44 
45 // DQM
46 #include "DQMServices/Core/interface/DQMEDAnalyzer.h"
47 #include "DQMServices/Core/interface/DQMStore.h"
48 
49 #include "DataFormats/Math/interface/deltaR.h"
50 #include "DataFormats/Math/interface/GeantUnits.h"
51 #include "CLHEP/Units/PhysicalConstants.h"
52 
53 //class declaration
54 class Primary4DVertexValidation : public DQMEDAnalyzer {
55  typedef math::XYZTLorentzVector LorentzVector;
56 
57  // auxiliary class holding simulated vertices
58  struct simPrimaryVertex {
59  simPrimaryVertex(double x1, double y1, double z1, double t1)
60  : x(x1),
61  y(y1),
62  z(z1),
63  t(t1),
64  ptsq(0),
65  closest_vertex_distance_z(-1.),
66  nGenTrk(0),
67  num_matched_reco_tracks(0),
68  average_match_quality(0.0) {
69  ptot.setPx(0);
70  ptot.setPy(0);
71  ptot.setPz(0);
72  ptot.setE(0);
73  p4 = LorentzVector(0, 0, 0, 0);
74  r = sqrt(x * x + y * y);
75  };
76  double x, y, z, r, t;
77  HepMC::FourVector ptot;
78  LorentzVector p4;
79  double ptsq;
80  double closest_vertex_distance_z;
81  int nGenTrk;
82  int num_matched_reco_tracks;
83  float average_match_quality;
84  EncodedEventId eventId;
85  TrackingVertexRef sim_vertex;
86  int OriginalIndex = -1;
87 
88  unsigned int nwosmatch = 0; // number of recvertices dominated by this simevt (by wos)
89  unsigned int nwntmatch = 0; // number of recvertices dominated by this simevt (by nt)
90  std::vector<unsigned int> wos_dominated_recv; // list of dominated recv (by wos, size==nwosmatch)
91 
92  std::map<unsigned int, double> wnt; // weighted number of tracks in recvtx (by index)
93  std::map<unsigned int, double> wos; // sum of wos in recvtx (by index) // oops -> this was int before 04-22
94  double sumwos = 0; // sum of wos in any recvtx
95  double sumwnt = 0; // sum of weighted tracks
96  unsigned int rec = NOT_MATCHED; // best match (NO_MATCH if not matched)
97  unsigned int matchQuality = 0; // quality flag
98 
99  void addTrack(unsigned int irecv, double twos, double twt) {
100  sumwnt += twt;
101  if (wnt.find(irecv) == wnt.end()) {
102  wnt[irecv] = twt;
103  } else {
104  wnt[irecv] += twt;
105  }
106 
107  sumwos += twos;
108  if (wos.find(irecv) == wos.end()) {
109  wos[irecv] = twos;
110  } else {
111  wos[irecv] += twos;
112  }
113  }
114  };
115 
116  // auxiliary class holding reconstructed vertices
117  struct recoPrimaryVertex {
118  recoPrimaryVertex(double x1, double y1, double z1)
119  : x(x1),
120  y(y1),
121  z(z1),
122  pt(0),
123  ptsq(0),
124  closest_vertex_distance_z(-1.),
125  nRecoTrk(0),
126  num_matched_sim_tracks(0),
127  ndof(0.),
128  recVtx(nullptr) {
129  r = sqrt(x * x + y * y);
130  };
131  double x, y, z, r;
132  double pt;
133  double ptsq;
134  double closest_vertex_distance_z;
135  int nRecoTrk;
136  int num_matched_sim_tracks;
137  double ndof;
138  const reco::Vertex* recVtx;
139  reco::VertexBaseRef recVtxRef;
140  int OriginalIndex = -1;
141 
142  std::map<unsigned int, double> wos; // simevent -> wos
143  std::map<unsigned int, double> wnt; // simevent -> weighted number of truth matched tracks
144  unsigned int wosmatch; // index of the simevent providing the largest contribution to wos
145  unsigned int wntmatch; // index of the simevent providing the highest number of tracks
146  double sumwos = 0; // total sum of wos of all truth matched tracks
147  double sumwnt = 0; // total weighted number of truth matchted tracks
148  double maxwos = 0; // largest wos sum from one sim event (wosmatch)
149  double maxwnt = 0; // largest weighted number of tracks from one sim event (ntmatch)
150  int maxwosnt = 0; // number of tracks from the simevt with highest wos
151  unsigned int sim = NOT_MATCHED; // best match (NO_MATCH if not matched)
152  unsigned int matchQuality = 0; // quality flag
153 
154  bool is_real() { return (matchQuality > 0) && (matchQuality < 99); }
155 
156  bool is_fake() { return (matchQuality <= 0) || (matchQuality >= 99); }
157 
158  bool is_signal() { return (sim == 0); }
159 
160  int split_from() {
161  if (is_real())
162  return -1;
163  if ((maxwos > 0) && (maxwos > 0.3 * sumwos))
164  return wosmatch;
165  return -1;
166  }
167  bool other_fake() { return (is_fake() && (split_from() < 0)); }
168 
169  void addTrack(unsigned int iev, double twos, double wt) {
170  sumwnt += wt;
171  if (wnt.find(iev) == wnt.end()) {
172  wnt[iev] = wt;
173  } else {
174  wnt[iev] += wt;
175  }
176 
177  sumwos += twos;
178  if (wos.find(iev) == wos.end()) {
179  wos[iev] = twos;
180  } else {
181  wos[iev] += twos;
182  }
183  }
184  };
185 
186 public:
187  explicit Primary4DVertexValidation(const edm::ParameterSet&);
188  ~Primary4DVertexValidation() override;
189 
190  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
191  void analyze(const edm::Event&, const edm::EventSetup&) override;
192  void bookHistograms(DQMStore::IBooker& i, edm::Run const&, edm::EventSetup const&) override;
193 
194 private:
195  void matchReco2Sim(std::vector<recoPrimaryVertex>&,
196  std::vector<simPrimaryVertex>&,
197  const edm::ValueMap<float>&,
198  const edm::ValueMap<float>&,
199  const edm::Handle<reco::BeamSpot>&);
200  bool matchRecoTrack2SimSignal(const reco::TrackBaseRef&);
201  const edm::Ref<std::vector<TrackingParticle>>* getMatchedTP(const reco::TrackBaseRef&, const TrackingVertexRef&);
202  double timeFromTrueMass(double, double, double, double);
203  bool select(const reco::Vertex&, int level = 0);
204  std::vector<Primary4DVertexValidation::simPrimaryVertex> getSimPVs(const edm::Handle<TrackingVertexCollection>&);
205  std::vector<Primary4DVertexValidation::recoPrimaryVertex> getRecoPVs(const edm::Handle<edm::View<reco::Vertex>>&);
206 
207  const bool mvaTPSel(const TrackingParticle&);
208  const bool mvaRecSel(const reco::TrackBase&, const reco::Vertex&, const double&, const double&);
209 
210  // ----------member data ---------------------------
211 
212  const std::string folder_;
213  static constexpr unsigned int NOT_MATCHED = 66666;
214  static constexpr double simUnit_ = 1e9; //sim time in s while reco time in ns
215  static constexpr double c_ = 2.99792458e1; //c in cm/ns
216  static constexpr double mvaL_ = 0.5; //MVA cuts for MVA categories
217  static constexpr double mvaH_ = 0.8;
218  static constexpr double selNdof_ = 4.;
219  static constexpr double maxRank_ = 8.;
220  static constexpr double maxTry_ = 10.;
221  static constexpr double zWosMatchMax_ = 1.;
222  static constexpr double etacutGEN_ = 4.; // |eta| < 4;
223  static constexpr double etacutREC_ = 3.; // |eta| < 3;
224  static constexpr double pTcut_ = 0.7; // PT > 0.7 GeV
225  static constexpr double deltaZcut_ = 0.1; // dz separation 1 mm
226  static constexpr double trackMaxBtlEta_ = 1.5;
227  static constexpr double trackMinEtlEta_ = 1.6;
228  static constexpr double trackMaxEtlEta_ = 3.;
229  static constexpr double tol_ = 1.e-4; // tolerance on reconstructed track time, [ns]
230 
231  static constexpr float c_cm_ns = geant_units::operators::convertMmToCm(CLHEP::c_light); // [mm/ns] -> [cm/ns]
232 
233  bool use_only_charged_tracks_;
234  bool debug_;
235  bool optionalPlots_;
236 
237  const double minProbHeavy_;
238  const double trackweightTh_;
239  const double mvaTh_;
240  const std::vector<double> lineDensityPar_;
241  const reco::RecoToSimCollection* r2s_;
242  const reco::SimToRecoCollection* s2r_;
243 
244  edm::EDGetTokenT<reco::TrackCollection> RecTrackToken_;
245 
246  edm::EDGetTokenT<std::vector<PileupSummaryInfo>> vecPileupSummaryInfoToken_;
247 
248  edm::EDGetTokenT<TrackingParticleCollection> trackingParticleCollectionToken_;
249  edm::EDGetTokenT<TrackingVertexCollection> trackingVertexCollectionToken_;
250  edm::EDGetTokenT<reco::SimToRecoCollection> simToRecoAssociationToken_;
251  edm::EDGetTokenT<reco::RecoToSimCollection> recoToSimAssociationToken_;
252  edm::EDGetTokenT<reco::BeamSpot> RecBeamSpotToken_;
253  edm::EDGetTokenT<edm::View<reco::Vertex>> Rec4DVerToken_;
254 
255  edm::EDGetTokenT<edm::ValueMap<int>> trackAssocToken_;
256  edm::EDGetTokenT<edm::ValueMap<float>> pathLengthToken_;
257  edm::EDGetTokenT<edm::ValueMap<float>> momentumToken_;
258  edm::EDGetTokenT<edm::ValueMap<float>> timeToken_;
259 
260  edm::EDGetTokenT<edm::ValueMap<float>> t0PidToken_;
261  edm::EDGetTokenT<edm::ValueMap<float>> t0SafePidToken_;
262  edm::EDGetTokenT<edm::ValueMap<float>> sigmat0SafePidToken_;
263  edm::EDGetTokenT<edm::ValueMap<float>> trackMVAQualToken_;
264 
265  edm::EDGetTokenT<edm::ValueMap<float>> tmtdToken_;
266  edm::EDGetTokenT<edm::ValueMap<float>> tofPiToken_;
267  edm::EDGetTokenT<edm::ValueMap<float>> tofKToken_;
268  edm::EDGetTokenT<edm::ValueMap<float>> tofPToken_;
269  edm::EDGetTokenT<edm::ValueMap<float>> probPiToken_;
270  edm::EDGetTokenT<edm::ValueMap<float>> probKToken_;
271  edm::EDGetTokenT<edm::ValueMap<float>> probPToken_;
272 
273  //histogram declaration
274  MonitorElement* meMVATrackEffPtTot_;
275  MonitorElement* meMVATrackMatchedEffPtTot_;
276  MonitorElement* meMVATrackMatchedEffPtMtd_;
277  MonitorElement* meMVATrackEffEtaTot_;
278  MonitorElement* meMVATrackMatchedEffEtaTot_;
279  MonitorElement* meMVATrackMatchedEffEtaMtd_;
280  MonitorElement* meMVATrackResTot_;
281  MonitorElement* meMVATrackPullTot_;
282  MonitorElement* meTrackResTot_;
283  MonitorElement* meTrackPullTot_;
284  MonitorElement* meTrackRes_[3];
285  MonitorElement* meTrackPull_[3];
286  MonitorElement* meTrackResMass_[3];
287  MonitorElement* meTrackResMassTrue_[3];
288  MonitorElement* meMVATrackZposResTot_;
289  MonitorElement* meTrackZposResTot_;
290  MonitorElement* meTrackZposRes_[3];
291  MonitorElement* meTrack3DposRes_[3];
292  MonitorElement* meTimeRes_;
293  MonitorElement* meTimePull_;
294  MonitorElement* meTimeSignalRes_;
295  MonitorElement* meTimeSignalPull_;
296  MonitorElement* mePUvsRealV_;
297  MonitorElement* mePUvsOtherFakeV_;
298  MonitorElement* mePUvsSplitV_;
299  MonitorElement* meMatchQual_;
300  MonitorElement* meDeltaZrealreal_;
301  MonitorElement* meDeltaZfakefake_;
302  MonitorElement* meDeltaZfakereal_;
303  MonitorElement* meDeltaTrealreal_;
304  MonitorElement* meDeltaTfakefake_;
305  MonitorElement* meDeltaTfakereal_;
306  MonitorElement* meRecoPosInSimCollection_;
307  MonitorElement* meRecoPosInRecoOrigCollection_;
308  MonitorElement* meSimPosInSimOrigCollection_;
309  MonitorElement* meRecoPVPosSignal_;
310  MonitorElement* meRecoPVPosSignalNotHighestPt_;
311  MonitorElement* meRecoVtxVsLineDensity_;
312  MonitorElement* meRecVerNumber_;
313  MonitorElement* meRecPVZ_;
314  MonitorElement* meRecPVT_;
315  MonitorElement* meSimPVZ_;
316 
317  //some tests
318  MonitorElement* meTrackResLowPTot_;
319  MonitorElement* meTrackResHighPTot_;
320  MonitorElement* meTrackPullLowPTot_;
321  MonitorElement* meTrackPullHighPTot_;
322 
323  MonitorElement* meTrackResLowP_[3];
324  MonitorElement* meTrackResHighP_[3];
325  MonitorElement* meTrackPullLowP_[3];
326  MonitorElement* meTrackPullHighP_[3];
327 
328  MonitorElement* meTrackResMassProtons_[3];
329  MonitorElement* meTrackResMassTrueProtons_[3];
330  MonitorElement* meTrackResMassPions_[3];
331  MonitorElement* meTrackResMassTruePions_[3];
332 
333  MonitorElement* meBarrelPIDp_;
334  MonitorElement* meEndcapPIDp_;
335 
336  MonitorElement* meBarrelNoPIDtype_;
337  MonitorElement* meEndcapNoPIDtype_;
338 
339  MonitorElement* meBarrelTruePiNoPID_;
340  MonitorElement* meBarrelTrueKNoPID_;
341  MonitorElement* meBarrelTruePNoPID_;
342  MonitorElement* meEndcapTruePiNoPID_;
343  MonitorElement* meEndcapTrueKNoPID_;
344  MonitorElement* meEndcapTruePNoPID_;
345 
346  MonitorElement* meBarrelTruePiAsPi_;
347  MonitorElement* meBarrelTruePiAsK_;
348  MonitorElement* meBarrelTruePiAsP_;
349  MonitorElement* meEndcapTruePiAsPi_;
350  MonitorElement* meEndcapTruePiAsK_;
351  MonitorElement* meEndcapTruePiAsP_;
352 
353  MonitorElement* meBarrelTrueKAsPi_;
354  MonitorElement* meBarrelTrueKAsK_;
355  MonitorElement* meBarrelTrueKAsP_;
356  MonitorElement* meEndcapTrueKAsPi_;
357  MonitorElement* meEndcapTrueKAsK_;
358  MonitorElement* meEndcapTrueKAsP_;
359 
360  MonitorElement* meBarrelTruePAsPi_;
361  MonitorElement* meBarrelTruePAsK_;
362  MonitorElement* meBarrelTruePAsP_;
363  MonitorElement* meEndcapTruePAsPi_;
364  MonitorElement* meEndcapTruePAsK_;
365  MonitorElement* meEndcapTruePAsP_;
366 };
367 
368 // constructors and destructor
369 Primary4DVertexValidation::Primary4DVertexValidation(const edm::ParameterSet& iConfig)
370  : folder_(iConfig.getParameter<std::string>("folder")),
371  use_only_charged_tracks_(iConfig.getParameter<bool>("useOnlyChargedTracks")),
372  debug_(iConfig.getUntrackedParameter<bool>("debug")),
373  optionalPlots_(iConfig.getUntrackedParameter<bool>("optionalPlots")),
374  minProbHeavy_(iConfig.getParameter<double>("minProbHeavy")),
375  trackweightTh_(iConfig.getParameter<double>("trackweightTh")),
376  mvaTh_(iConfig.getParameter<double>("mvaTh")),
377  lineDensityPar_(iConfig.getParameter<std::vector<double>>("lineDensityPar")) {
378  vecPileupSummaryInfoToken_ = consumes<std::vector<PileupSummaryInfo>>(edm::InputTag(std::string("addPileupInfo")));
379  trackingParticleCollectionToken_ =
380  consumes<TrackingParticleCollection>(iConfig.getParameter<edm::InputTag>("SimTag"));
381  trackingVertexCollectionToken_ = consumes<TrackingVertexCollection>(iConfig.getParameter<edm::InputTag>("SimTag"));
382  simToRecoAssociationToken_ =
383  consumes<reco::SimToRecoCollection>(iConfig.getParameter<edm::InputTag>("TPtoRecoTrackAssoc"));
384  recoToSimAssociationToken_ =
385  consumes<reco::RecoToSimCollection>(iConfig.getParameter<edm::InputTag>("TPtoRecoTrackAssoc"));
386  RecTrackToken_ = consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("mtdTracks"));
387  RecBeamSpotToken_ = consumes<reco::BeamSpot>(iConfig.getParameter<edm::InputTag>("offlineBS"));
388  Rec4DVerToken_ = consumes<edm::View<reco::Vertex>>(iConfig.getParameter<edm::InputTag>("offline4DPV"));
389  trackAssocToken_ = consumes<edm::ValueMap<int>>(iConfig.getParameter<edm::InputTag>("trackAssocSrc"));
390  pathLengthToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("pathLengthSrc"));
391  momentumToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("momentumSrc"));
392  timeToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("timeSrc"));
393  t0PidToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("t0PID"));
394  t0SafePidToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("t0SafePID"));
395  sigmat0SafePidToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("sigmat0SafePID"));
396  trackMVAQualToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("trackMVAQual"));
397  tmtdToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("tmtd"));
398  tofPiToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("tofPi"));
399  tofKToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("tofK"));
400  tofPToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("tofP"));
401  probPiToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("probPi"));
402  probKToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("probK"));
403  probPToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("probP"));
404 }
405 
406 Primary4DVertexValidation::~Primary4DVertexValidation() {}
407 
408 //
409 // member functions
410 //
411 void Primary4DVertexValidation::bookHistograms(DQMStore::IBooker& ibook,
412  edm::Run const& iRun,
413  edm::EventSetup const& iSetup) {
414  ibook.setCurrentFolder(folder_);
415  // --- histograms booking
416  meMVATrackEffPtTot_ = ibook.book1D("MVAEffPtTot", "Pt of tracks associated to LV; track pt [GeV] ", 110, 0., 11.);
417  meMVATrackEffEtaTot_ = ibook.book1D("MVAEffEtaTot", "Pt of tracks associated to LV; track eta ", 66, 0., 3.3);
418  meMVATrackMatchedEffPtTot_ =
419  ibook.book1D("MVAMatchedEffPtTot", "Pt of tracks associated to LV matched to TP; track pt [GeV] ", 110, 0., 11.);
420  meMVATrackMatchedEffPtMtd_ = ibook.book1D(
421  "MVAMatchedEffPtMtd", "Pt of tracks associated to LV matched to TP with time; track pt [GeV] ", 110, 0., 11.);
422  meMVATrackMatchedEffEtaTot_ =
423  ibook.book1D("MVAMatchedEffEtaTot", "Pt of tracks associated to LV matched to TP; track eta ", 66, 0., 3.3);
424  meMVATrackMatchedEffEtaMtd_ = ibook.book1D(
425  "MVAMatchedEffEtaMtd", "Pt of tracks associated to LV matched to TP with time; track eta ", 66, 0., 3.3);
426  meMVATrackResTot_ = ibook.book1D(
427  "MVATrackRes", "t_{rec} - t_{sim} for tracks from LV MVA sel.; t_{rec} - t_{sim} [ns] ", 120, -0.15, 0.15);
428  meTrackResTot_ = ibook.book1D("TrackRes", "t_{rec} - t_{sim} for tracks; t_{rec} - t_{sim} [ns] ", 120, -0.15, 0.15);
429  meTrackRes_[0] = ibook.book1D(
430  "TrackRes-LowMVA", "t_{rec} - t_{sim} for tracks with MVA < 0.5; t_{rec} - t_{sim} [ns] ", 100, -1., 1.);
431  meTrackRes_[1] = ibook.book1D(
432  "TrackRes-MediumMVA", "t_{rec} - t_{sim} for tracks with 0.5 < MVA < 0.8; t_{rec} - t_{sim} [ns] ", 100, -1., 1.);
433  meTrackRes_[2] = ibook.book1D(
434  "TrackRes-HighMVA", "t_{rec} - t_{sim} for tracks with MVA > 0.8; t_{rec} - t_{sim} [ns] ", 100, -1., 1.);
435  if (optionalPlots_) {
436  meTrackResMass_[0] = ibook.book1D(
437  "TrackResMass-LowMVA", "t_{rec} - t_{est} for tracks with MVA < 0.5; t_{rec} - t_{est} [ns] ", 100, -1., 1.);
438  meTrackResMass_[1] = ibook.book1D("TrackResMass-MediumMVA",
439  "t_{rec} - t_{est} for tracks with 0.5 < MVA < 0.8; t_{rec} - t_{est} [ns] ",
440  100,
441  -1.,
442  1.);
443  meTrackResMass_[2] = ibook.book1D(
444  "TrackResMass-HighMVA", "t_{rec} - t_{est} for tracks with MVA > 0.8; t_{rec} - t_{est} [ns] ", 100, -1., 1.);
445  meTrackResMassTrue_[0] = ibook.book1D(
446  "TrackResMassTrue-LowMVA", "t_{est} - t_{sim} for tracks with MVA < 0.5; t_{est} - t_{sim} [ns] ", 100, -1., 1.);
447  meTrackResMassTrue_[1] = ibook.book1D("TrackResMassTrue-MediumMVA",
448  "t_{est} - t_{sim} for tracks with 0.5 < MVA < 0.8; t_{est} - t_{sim} [ns] ",
449  100,
450  -1.,
451  1.);
452  meTrackResMassTrue_[2] = ibook.book1D("TrackResMassTrue-HighMVA",
453  "t_{est} - t_{sim} for tracks with MVA > 0.8; t_{est} - t_{sim} [ns] ",
454  100,
455  -1.,
456  1.);
457  }
458  meMVATrackPullTot_ =
459  ibook.book1D("MVATrackPull", "Pull for tracks from LV MAV sel.; (t_{rec}-t_{sim})/#sigma_{t}", 50, -5., 5.);
460  meTrackPullTot_ = ibook.book1D("TrackPull", "Pull for tracks; (t_{rec}-t_{sim})/#sigma_{t}", 100, -10., 10.);
461  meTrackPull_[0] =
462  ibook.book1D("TrackPull-LowMVA", "Pull for tracks with MVA < 0.5; (t_{rec}-t_{sim})/#sigma_{t}", 100, -10., 10.);
463  meTrackPull_[1] = ibook.book1D(
464  "TrackPull-MediumMVA", "Pull for tracks with 0.5 < MVA < 0.8; (t_{rec}-t_{sim})/#sigma_{t}", 100, -10., 10.);
465  meTrackPull_[2] =
466  ibook.book1D("TrackPull-HighMVA", "Pull for tracks with MVA > 0.8; (t_{rec}-t_{sim})/#sigma_{t}", 100, -10., 10.);
467  meMVATrackZposResTot_ = ibook.book1D(
468  "MVATrackZposResTot", "Z_{PCA} - Z_{sim} for tracks from LV MVA sel.;Z_{PCA} - Z_{sim} [cm] ", 50, -0.1, 0.1);
469  meTrackZposResTot_ =
470  ibook.book1D("TrackZposResTot", "Z_{PCA} - Z_{sim} for tracks;Z_{PCA} - Z_{sim} [cm] ", 50, -0.5, 0.5);
471  meTrackZposRes_[0] = ibook.book1D(
472  "TrackZposRes-LowMVA", "Z_{PCA} - Z_{sim} for tracks with MVA < 0.5;Z_{PCA} - Z_{sim} [cm] ", 50, -0.5, 0.5);
473  meTrackZposRes_[1] = ibook.book1D("TrackZposRes-MediumMVA",
474  "Z_{PCA} - Z_{sim} for tracks with 0.5 < MVA < 0.8 ;Z_{PCA} - Z_{sim} [cm] ",
475  50,
476  -0.5,
477  0.5);
478  meTrackZposRes_[2] = ibook.book1D(
479  "TrackZposRes-HighMVA", "Z_{PCA} - Z_{sim} for tracks with MVA > 0.8 ;Z_{PCA} - Z_{sim} [cm] ", 50, -0.5, 0.5);
480  meTrack3DposRes_[0] =
481  ibook.book1D("Track3DposRes-LowMVA",
482  "3dPos_{PCA} - 3dPos_{sim} for tracks with MVA < 0.5 ;3dPos_{PCA} - 3dPos_{sim} [cm] ",
483  50,
484  -0.5,
485  0.5);
486  meTrack3DposRes_[1] =
487  ibook.book1D("Track3DposRes-MediumMVA",
488  "3dPos_{PCA} - 3dPos_{sim} for tracks with 0.5 < MVA < 0.8 ;3dPos_{PCA} - 3dPos_{sim} [cm] ",
489  50,
490  -0.5,
491  0.5);
492  meTrack3DposRes_[2] =
493  ibook.book1D("Track3DposRes-HighMVA",
494  "3dPos_{PCA} - 3dPos_{sim} for tracks with MVA > 0.8;3dPos_{PCA} - 3dPos_{sim} [cm] ",
495  50,
496  -0.5,
497  0.5);
498  meTimeRes_ = ibook.book1D("TimeRes", "t_{rec} - t_{sim} ;t_{rec} - t_{sim} [ns] ", 40, -0.2, 0.2);
499  meTimePull_ = ibook.book1D("TimePull", "Pull; t_{rec} - t_{sim}/#sigma_{t rec}", 100, -10., 10.);
500  meTimeSignalRes_ =
501  ibook.book1D("TimeSignalRes", "t_{rec} - t_{sim} for signal ;t_{rec} - t_{sim} [ns] ", 40, -0.2, 0.2);
502  meTimeSignalPull_ =
503  ibook.book1D("TimeSignalPull", "Pull for signal; t_{rec} - t_{sim}/#sigma_{t rec}", 100, -10., 10.);
504  mePUvsRealV_ =
505  ibook.bookProfile("PUvsReal", "#PU vertices vs #real matched vertices;#PU;#real ", 100, 0, 300, 100, 0, 200);
506  mePUvsOtherFakeV_ = ibook.bookProfile(
507  "PUvsOtherFake", "#PU vertices vs #other fake matched vertices;#PU;#other fake ", 100, 0, 300, 100, 0, 20);
508  mePUvsSplitV_ =
509  ibook.bookProfile("PUvsSplit", "#PU vertices vs #split matched vertices;#PU;#split ", 100, 0, 300, 100, 0, 20);
510  meMatchQual_ = ibook.book1D("MatchQuality", "RECO-SIM vertex match quality; ", 8, 0, 8.);
511  meDeltaZrealreal_ = ibook.book1D("DeltaZrealreal", "#Delta Z real-real; |#Delta Z (r-r)| [cm]", 100, 0, 0.5);
512  meDeltaZfakefake_ = ibook.book1D("DeltaZfakefake", "#Delta Z fake-fake; |#Delta Z (f-f)| [cm]", 100, 0, 0.5);
513  meDeltaZfakereal_ = ibook.book1D("DeltaZfakereal", "#Delta Z fake-real; |#Delta Z (f-r)| [cm]", 100, 0, 0.5);
514  meDeltaTrealreal_ = ibook.book1D("DeltaTrealreal", "#Delta T real-real; |#Delta T (r-r)| [sigma]", 60, 0., 30.);
515  meDeltaTfakefake_ = ibook.book1D("DeltaTfakefake", "#Delta T fake-fake; |#Delta T (f-f)| [sigma]", 60, 0., 30.);
516  meDeltaTfakereal_ = ibook.book1D("DeltaTfakereal", "#Delta T fake-real; |#Delta T (f-r)| [sigma]", 60, 0., 30.);
517  if (optionalPlots_) {
518  meRecoPosInSimCollection_ = ibook.book1D(
519  "RecoPosInSimCollection", "Sim signal vertex index associated to Reco signal vertex; Sim PV index", 200, 0, 200);
520  meRecoPosInRecoOrigCollection_ =
521  ibook.book1D("RecoPosInRecoOrigCollection", "Reco signal index in OrigCollection; Reco index", 200, 0, 200);
522  meSimPosInSimOrigCollection_ =
523  ibook.book1D("SimPosInSimOrigCollection", "Sim signal index in OrigCollection; Sim index", 200, 0, 200);
524  }
525  meRecoPVPosSignal_ =
526  ibook.book1D("RecoPVPosSignal", "Position in reco collection of PV associated to sim signal", 200, 0, 200);
527  meRecoPVPosSignalNotHighestPt_ =
528  ibook.book1D("RecoPVPosSignalNotHighestPt",
529  "Position in reco collection of PV associated to sim signal not highest Pt",
530  200,
531  0,
532  200);
533  meRecoVtxVsLineDensity_ =
534  ibook.book1D("RecoVtxVsLineDensity", "#Reco vertices/mm/event; line density [#vtx/mm/event]", 160, 0., 4.);
535  meRecVerNumber_ = ibook.book1D("RecVerNumber", "RECO Vertex Number: Number of vertices", 50, 0, 250);
536  meRecPVZ_ = ibook.book1D("recPVZ", "Weighted #Rec vertices/mm", 400, -20., 20.);
537  meRecPVT_ = ibook.book1D("recPVT", "#Rec vertices/10 ps", 200, -1., 1.);
538  meSimPVZ_ = ibook.book1D("simPVZ", "Weighted #Sim vertices/mm", 400, -20., 20.);
539 
540  //some tests
541  meTrackResLowPTot_ = ibook.book1D(
542  "TrackResLowP", "t_{rec} - t_{sim} for tracks with p < 2 GeV; t_{rec} - t_{sim} [ns] ", 70, -0.15, 0.15);
543  meTrackResLowP_[0] =
544  ibook.book1D("TrackResLowP-LowMVA",
545  "t_{rec} - t_{sim} for tracks with MVA < 0.5 and p < 2 GeV; t_{rec} - t_{sim} [ns] ",
546  100,
547  -1.,
548  1.);
549  meTrackResLowP_[1] =
550  ibook.book1D("TrackResLowP-MediumMVA",
551  "t_{rec} - t_{sim} for tracks with 0.5 < MVA < 0.8 and p < 2 GeV; t_{rec} - t_{sim} [ns] ",
552  100,
553  -1.,
554  1.);
555  meTrackResLowP_[2] =
556  ibook.book1D("TrackResLowP-HighMVA",
557  "t_{rec} - t_{sim} for tracks with MVA > 0.8 and p < 2 GeV; t_{rec} - t_{sim} [ns] ",
558  100,
559  -1.,
560  1.);
561  meTrackResHighPTot_ = ibook.book1D(
562  "TrackResHighP", "t_{rec} - t_{sim} for tracks with p > 2 GeV; t_{rec} - t_{sim} [ns] ", 70, -0.15, 0.15);
563  meTrackResHighP_[0] =
564  ibook.book1D("TrackResHighP-LowMVA",
565  "t_{rec} - t_{sim} for tracks with MVA < 0.5 and p > 2 GeV; t_{rec} - t_{sim} [ns] ",
566  100,
567  -1.,
568  1.);
569  meTrackResHighP_[1] =
570  ibook.book1D("TrackResHighP-MediumMVA",
571  "t_{rec} - t_{sim} for tracks with 0.5 < MVA < 0.8 and p > 2 GeV; t_{rec} - t_{sim} [ns] ",
572  100,
573  -1.,
574  1.);
575  meTrackResHighP_[2] =
576  ibook.book1D("TrackResHighP-HighMVA",
577  "t_{rec} - t_{sim} for tracks with MVA > 0.8 and p > 2 GeV; t_{rec} - t_{sim} [ns] ",
578  100,
579  -1.,
580  1.);
581  meTrackPullLowPTot_ =
582  ibook.book1D("TrackPullLowP", "Pull for tracks with p < 2 GeV; (t_{rec}-t_{sim})/#sigma_{t}", 100, -10., 10.);
583  meTrackPullLowP_[0] = ibook.book1D("TrackPullLowP-LowMVA",
584  "Pull for tracks with MVA < 0.5 and p < 2 GeV; (t_{rec}-t_{sim})/#sigma_{t}",
585  100,
586  -10.,
587  10.);
588  meTrackPullLowP_[1] = ibook.book1D("TrackPullLowP-MediumMVA",
589  "Pull for tracks with 0.5 < MVA < 0.8 and p < 2 GeV; (t_{rec}-t_{sim})/#sigma_{t}",
590  100,
591  -10.,
592  10.);
593  meTrackPullLowP_[2] = ibook.book1D("TrackPullLowP-HighMVA",
594  "Pull for tracks with MVA > 0.8 and p < 2 GeV; (t_{rec}-t_{sim})/#sigma_{t}",
595  100,
596  -10.,
597  10.);
598  meTrackPullHighPTot_ =
599  ibook.book1D("TrackPullHighP", "Pull for tracks with p > 2 GeV; (t_{rec}-t_{sim})/#sigma_{t}", 100, -10., 10.);
600  meTrackPullHighP_[0] = ibook.book1D("TrackPullHighP-LowMVA",
601  "Pull for tracks with MVA < 0.5 and p > 2 GeV; (t_{rec}-t_{sim})/#sigma_{t}",
602  100,
603  -10.,
604  10.);
605  meTrackPullHighP_[1] =
606  ibook.book1D("TrackPullHighP-MediumMVA",
607  "Pull for tracks with 0.5 < MVA < 0.8 and p > 2 GeV; (t_{rec}-t_{sim})/#sigma_{t}",
608  100,
609  -10.,
610  10.);
611  meTrackPullHighP_[2] = ibook.book1D("TrackPullHighP-HighMVA",
612  "Pull for tracks with MVA > 0.8 and p > 2 GeV; (t_{rec}-t_{sim})/#sigma_{t}",
613  100,
614  -10.,
615  10.);
616  if (optionalPlots_) {
617  meTrackResMassProtons_[0] =
618  ibook.book1D("TrackResMass-Protons-LowMVA",
619  "t_{rec} - t_{est} for proton tracks with MVA < 0.5; t_{rec} - t_{est} [ns] ",
620  100,
621  -1.,
622  1.);
623  meTrackResMassProtons_[1] =
624  ibook.book1D("TrackResMass-Protons-MediumMVA",
625  "t_{rec} - t_{est} for proton tracks with 0.5 < MVA < 0.8; t_{rec} - t_{est} [ns] ",
626  100,
627  -1.,
628  1.);
629  meTrackResMassProtons_[2] =
630  ibook.book1D("TrackResMass-Protons-HighMVA",
631  "t_{rec} - t_{est} for proton tracks with MVA > 0.8; t_{rec} - t_{est} [ns] ",
632  100,
633  -1.,
634  1.);
635  meTrackResMassTrueProtons_[0] =
636  ibook.book1D("TrackResMassTrue-Protons-LowMVA",
637  "t_{est} - t_{sim} for proton tracks with MVA < 0.5; t_{est} - t_{sim} [ns] ",
638  100,
639  -1.,
640  1.);
641  meTrackResMassTrueProtons_[1] =
642  ibook.book1D("TrackResMassTrue-Protons-MediumMVA",
643  "t_{est} - t_{sim} for proton tracks with 0.5 < MVA < 0.8; t_{est} - t_{sim} [ns] ",
644  100,
645  -1.,
646  1.);
647  meTrackResMassTrueProtons_[2] =
648  ibook.book1D("TrackResMassTrue-Protons-HighMVA",
649  "t_{est} - t_{sim} for proton tracks with MVA > 0.8; t_{est} - t_{sim} [ns] ",
650  100,
651  -1.,
652  1.);
653 
654  meTrackResMassPions_[0] = ibook.book1D("TrackResMass-Pions-LowMVA",
655  "t_{rec} - t_{est} for pion tracks with MVA < 0.5; t_{rec} - t_{est} [ns] ",
656  100,
657  -1.,
658  1.);
659  meTrackResMassPions_[1] =
660  ibook.book1D("TrackResMass-Pions-MediumMVA",
661  "t_{rec} - t_{est} for pion tracks with 0.5 < MVA < 0.8; t_{rec} - t_{est} [ns] ",
662  100,
663  -1.,
664  1.);
665  meTrackResMassPions_[2] = ibook.book1D("TrackResMass-Pions-HighMVA",
666  "t_{rec} - t_{est} for pion tracks with MVA > 0.8; t_{rec} - t_{est} [ns] ",
667  100,
668  -1.,
669  1.);
670  meTrackResMassTruePions_[0] =
671  ibook.book1D("TrackResMassTrue-Pions-LowMVA",
672  "t_{est} - t_{sim} for pion tracks with MVA < 0.5; t_{est} - t_{sim} [ns] ",
673  100,
674  -1.,
675  1.);
676  meTrackResMassTruePions_[1] =
677  ibook.book1D("TrackResMassTrue-Pions-MediumMVA",
678  "t_{est} - t_{sim} for pion tracks with 0.5 < MVA < 0.8; t_{est} - t_{sim} [ns] ",
679  100,
680  -1.,
681  1.);
682  meTrackResMassTruePions_[2] =
683  ibook.book1D("TrackResMassTrue-Pions-HighMVA",
684  "t_{est} - t_{sim} for pion tracks with MVA > 0.8; t_{est} - t_{sim} [ns] ",
685  100,
686  -1.,
687  1.);
688  }
689 
690  meBarrelPIDp_ = ibook.book1D("BarrelPIDp", "PID track MTD momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
691  meEndcapPIDp_ = ibook.book1D("EndcapPIDp", "PID track with MTD momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
692 
693  meBarrelNoPIDtype_ = ibook.book1D("BarrelNoPIDtype", "Barrel PID failure category", 4, 0., 4.);
694  meEndcapNoPIDtype_ = ibook.book1D("EndcapNoPIDtype", "Endcap PID failure category", 4, 0., 4.);
695 
696  meBarrelNoPIDtype_ = ibook.book1D("BarrelNoPIDtype", "Barrel PID failure category", 4, 0., 4.);
697  meEndcapNoPIDtype_ = ibook.book1D("EndcapNoPIDtype", "Endcap PID failure category", 4, 0., 4.);
698 
699  meBarrelTruePiNoPID_ =
700  ibook.book1D("BarrelTruePiNoPID", "True pi NoPID momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
701  meBarrelTrueKNoPID_ =
702  ibook.book1D("BarrelTrueKNoPID", "True k NoPID momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
703  meBarrelTruePNoPID_ =
704  ibook.book1D("BarrelTruePNoPID", "True p NoPID momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
705  meEndcapTruePiNoPID_ =
706  ibook.book1D("EndcapTruePiNoPID", "True NoPIDpi momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
707  meEndcapTrueKNoPID_ =
708  ibook.book1D("EndcapTrueKNoPID", "True k NoPID momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
709  meEndcapTruePNoPID_ =
710  ibook.book1D("EndcapTruePNoPID", "True p NoPID momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
711 
712  meBarrelTruePiAsPi_ =
713  ibook.book1D("BarrelTruePiAsPi", "True pi as pi momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
714  meBarrelTruePiAsK_ =
715  ibook.book1D("BarrelTruePiAsK", "True pi as k momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
716  meBarrelTruePiAsP_ =
717  ibook.book1D("BarrelTruePiAsP", "True pi as p momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
718  meEndcapTruePiAsPi_ =
719  ibook.book1D("EndcapTruePiAsPi", "True pi as pi momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
720  meEndcapTruePiAsK_ =
721  ibook.book1D("EndcapTruePiAsK", "True pi as k momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
722  meEndcapTruePiAsP_ =
723  ibook.book1D("EndcapTruePiAsP", "True pi as p momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
724 
725  meBarrelTrueKAsPi_ =
726  ibook.book1D("BarrelTrueKAsPi", "True k as pi momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
727  meBarrelTrueKAsK_ = ibook.book1D("BarrelTrueKAsK", "True k as k momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
728  meBarrelTrueKAsP_ = ibook.book1D("BarrelTrueKAsP", "True k as p momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
729  meEndcapTrueKAsPi_ =
730  ibook.book1D("EndcapTrueKAsPi", "True k as pi momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
731  meEndcapTrueKAsK_ = ibook.book1D("EndcapTrueKAsK", "True k as k momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
732  meEndcapTrueKAsP_ = ibook.book1D("EndcapTrueKAsP", "True k as p momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
733 
734  meBarrelTruePAsPi_ =
735  ibook.book1D("BarrelTruePAsPi", "True p as pi momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
736  meBarrelTruePAsK_ = ibook.book1D("BarrelTruePAsK", "True p as k momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
737  meBarrelTruePAsP_ = ibook.book1D("BarrelTruePAsP", "True p as p momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
738  meEndcapTruePAsPi_ =
739  ibook.book1D("EndcapTruePAsPi", "True p as pi momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
740  meEndcapTruePAsK_ = ibook.book1D("EndcapTruePAsK", "True p as k momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
741  meEndcapTruePAsP_ = ibook.book1D("EndcapTruePAsP", "True p as p momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
742 }
743 
744 bool Primary4DVertexValidation::matchRecoTrack2SimSignal(const reco::TrackBaseRef& recoTrack) {
745  auto found = r2s_->find(recoTrack);
746 
747  // reco track not matched to any TP
748  if (found == r2s_->end())
749  return false;
750 
752  for (const auto& tp : found->val) {
753  if (tp.first->eventId().bunchCrossing() == 0 && tp.first->eventId().event() == 0)
754  return true;
755  }
756 
757  // reco track not matched to any TP from signal vertex
758  return false;
759 }
760 
761 const edm::Ref<std::vector<TrackingParticle>>* Primary4DVertexValidation::getMatchedTP(
762  const reco::TrackBaseRef& recoTrack, const TrackingVertexRef& vsim) {
763  auto found = r2s_->find(recoTrack);
764 
765  // reco track not matched to any TP
766  if (found == r2s_->end())
767  return nullptr;
768 
769  //matched TP equal to any TP of sim vertex
770  for (const auto& tp : found->val) {
771  if (std::find_if(vsim->daughterTracks_begin(), vsim->daughterTracks_end(), [&](const TrackingParticleRef& vtp) {
772  return tp.first == vtp;
773  }) != vsim->daughterTracks_end())
774  return &tp.first;
775  }
776 
777  // reco track not matched to any TP from vertex
778  return nullptr;
779 }
780 
781 double Primary4DVertexValidation::timeFromTrueMass(double mass, double pathlength, double momentum, double time) {
782  if (time > 0 && pathlength > 0 && mass > 0) {
783  double gammasq = 1. + momentum * momentum / (mass * mass);
784  double v = c_ * std::sqrt(1. - 1. / gammasq); // cm / ns
785  double t_est = time - (pathlength / v);
786 
787  return t_est;
788  } else {
789  return -1;
790  }
791 }
792 
793 bool Primary4DVertexValidation::select(const reco::Vertex& v, int level) {
794  /* level
795  0 !isFake && ndof>4 (default)
796  1 !isFake && ndof>4 && prob > 0.01
797  2 !isFake && ndof>4 && prob > 0.01 && ptmax2 > 0.4
798  */
799  if (v.isFake())
800  return false;
801  if ((level == 0) && (v.ndof() > selNdof_))
802  return true;
803  /*if ((level == 1) && (v.ndof() > selNdof_) && (vertex_pxy(v) > 0.01))
804  return true;
805  if ((level == 2) && (v.ndof() > selNdof_) && (vertex_pxy(v) > 0.01) && (vertex_ptmax2(v) > 0.4))
806  return true;
807  if ((level == 3) && (v.ndof() > selNdof_) && (vertex_ptmax2(v) < 0.4))
808  return true;*/
809  return false;
810 }
811 
812 /* Extract information form TrackingParticles/TrackingVertex and fill
813  * the helper class simPrimaryVertex with proper generation-level
814  * information */
815 std::vector<Primary4DVertexValidation::simPrimaryVertex> Primary4DVertexValidation::getSimPVs(
816  const edm::Handle<TrackingVertexCollection>& tVC) {
817  std::vector<Primary4DVertexValidation::simPrimaryVertex> simpv;
818  int current_event = -1;
819  int s = -1;
820  for (TrackingVertexCollection::const_iterator v = tVC->begin(); v != tVC->end(); ++v) {
821  //We keep only the first vertex from all the events at BX=0.
822  if (v->eventId().bunchCrossing() != 0)
823  continue;
824  if (v->eventId().event() != current_event) {
825  current_event = v->eventId().event();
826  } else {
827  continue;
828  }
829  s++;
830  if (std::abs(v->position().z()) > 1000)
831  continue; // skip junk vertices
832 
833  // could be a new vertex, check all primaries found so far to avoid multiple entries
834  simPrimaryVertex sv(v->position().x(), v->position().y(), v->position().z(), v->position().t());
835  sv.eventId = v->eventId();
836  sv.sim_vertex = TrackingVertexRef(tVC, std::distance(tVC->begin(), v));
837  sv.OriginalIndex = s;
838 
839  for (TrackingParticleRefVector::iterator iTrack = v->daughterTracks_begin(); iTrack != v->daughterTracks_end();
840  ++iTrack) {
841  assert((**iTrack).eventId().bunchCrossing() == 0);
842  }
843  simPrimaryVertex* vp = nullptr; // will become non-NULL if a vertex is found and then point to it
844  for (std::vector<simPrimaryVertex>::iterator v0 = simpv.begin(); v0 != simpv.end(); v0++) {
845  if ((sv.eventId == v0->eventId) && (std::abs(sv.x - v0->x) < 1e-5) && (std::abs(sv.y - v0->y) < 1e-5) &&
846  (std::abs(sv.z - v0->z) < 1e-5)) {
847  vp = &(*v0);
848  break;
849  }
850  }
851  if (!vp) {
852  // this is a new vertex, add it to the list of sim-vertices
853  simpv.push_back(sv);
854  vp = &simpv.back();
855  }
856 
857  // Loop over daughter track(s) as Tracking Particles
858  for (TrackingVertex::tp_iterator iTP = v->daughterTracks_begin(); iTP != v->daughterTracks_end(); ++iTP) {
859  auto momentum = (*(*iTP)).momentum();
860  const reco::Track* matched_best_reco_track = nullptr;
861  double match_quality = -1;
862  if (use_only_charged_tracks_ && (**iTP).charge() == 0)
863  continue;
864  if (s2r_->find(*iTP) != s2r_->end()) {
865  matched_best_reco_track = (*s2r_)[*iTP][0].first.get();
866  match_quality = (*s2r_)[*iTP][0].second;
867  }
868 
869  vp->ptot.setPx(vp->ptot.x() + momentum.x());
870  vp->ptot.setPy(vp->ptot.y() + momentum.y());
871  vp->ptot.setPz(vp->ptot.z() + momentum.z());
872  vp->ptot.setE(vp->ptot.e() + (**iTP).energy());
873  vp->ptsq += ((**iTP).pt() * (**iTP).pt());
874 
875  if (matched_best_reco_track) {
876  vp->num_matched_reco_tracks++;
877  vp->average_match_quality += match_quality;
878  }
879  } // End of for loop on daughters sim-particles
880  if (vp->num_matched_reco_tracks)
881  vp->average_match_quality /= static_cast<float>(vp->num_matched_reco_tracks);
882  if (debug_) {
883  edm::LogPrint("Primary4DVertexValidation")
884  << "average number of associated tracks: " << vp->num_matched_reco_tracks / static_cast<float>(vp->nGenTrk)
885  << " with average quality: " << vp->average_match_quality;
886  }
887  } // End of for loop on tracking vertices
888 
889  // In case of no simulated vertices, break here
890  if (simpv.empty())
891  return simpv;
892 
893  // Now compute the closest distance in z between all simulated vertex
894  // first initialize
895  auto prev_z = simpv.back().z;
896  for (simPrimaryVertex& vsim : simpv) {
897  vsim.closest_vertex_distance_z = std::abs(vsim.z - prev_z);
898  prev_z = vsim.z;
899  }
900  // then calculate
901  for (std::vector<simPrimaryVertex>::iterator vsim = simpv.begin(); vsim != simpv.end(); vsim++) {
902  std::vector<simPrimaryVertex>::iterator vsim2 = vsim;
903  vsim2++;
904  for (; vsim2 != simpv.end(); vsim2++) {
905  double distance = std::abs(vsim->z - vsim2->z);
906  // need both to be complete
907  vsim->closest_vertex_distance_z = std::min(vsim->closest_vertex_distance_z, distance);
908  vsim2->closest_vertex_distance_z = std::min(vsim2->closest_vertex_distance_z, distance);
909  }
910  }
911  return simpv;
912 }
913 
914 /* Extract information form recoVertex and fill the helper class
915  * recoPrimaryVertex with proper reco-level information */
916 std::vector<Primary4DVertexValidation::recoPrimaryVertex> Primary4DVertexValidation::getRecoPVs(
917  const edm::Handle<edm::View<reco::Vertex>>& tVC) {
918  std::vector<Primary4DVertexValidation::recoPrimaryVertex> recopv;
919  int r = -1;
920  for (auto v = tVC->begin(); v != tVC->end(); ++v) {
921  r++;
922  // Skip junk vertices
923  if (std::abs(v->z()) > 1000)
924  continue;
925  if (v->isFake() || !v->isValid())
926  continue;
927 
928  recoPrimaryVertex sv(v->position().x(), v->position().y(), v->position().z());
929  sv.recVtx = &(*v);
930  sv.recVtxRef = reco::VertexBaseRef(tVC, std::distance(tVC->begin(), v));
931 
932  sv.OriginalIndex = r;
933  sv.ndof = v->ndof();
934  // this is a new vertex, add it to the list of reco-vertices
935  recopv.push_back(sv);
936  Primary4DVertexValidation::recoPrimaryVertex* vp = &recopv.back();
937 
938  // Loop over daughter track(s)
939  for (auto iTrack = v->tracks_begin(); iTrack != v->tracks_end(); ++iTrack) {
940  auto momentum = (*(*iTrack)).innerMomentum();
941  if (momentum.mag2() == 0)
942  momentum = (*(*iTrack)).momentum();
943  vp->pt += std::sqrt(momentum.perp2());
944  vp->ptsq += (momentum.perp2());
945  vp->nRecoTrk++;
946 
947  auto matched = r2s_->find(*iTrack);
948  if (matched != r2s_->end()) {
949  vp->num_matched_sim_tracks++;
950  }
951 
952  } // End of for loop on daughters reconstructed tracks
953  } // End of for loop on tracking vertices
954 
955  // In case of no reco vertices, break here
956  if (recopv.empty())
957  return recopv;
958 
959  // Now compute the closest distance in z between all reconstructed vertex
960  // first initialize
961  auto prev_z = recopv.back().z;
962  for (recoPrimaryVertex& vreco : recopv) {
963  vreco.closest_vertex_distance_z = std::abs(vreco.z - prev_z);
964  prev_z = vreco.z;
965  }
966  for (std::vector<recoPrimaryVertex>::iterator vreco = recopv.begin(); vreco != recopv.end(); vreco++) {
967  std::vector<recoPrimaryVertex>::iterator vreco2 = vreco;
968  vreco2++;
969  for (; vreco2 != recopv.end(); vreco2++) {
970  double distance = std::abs(vreco->z - vreco2->z);
971  // need both to be complete
972  vreco->closest_vertex_distance_z = std::min(vreco->closest_vertex_distance_z, distance);
973  vreco2->closest_vertex_distance_z = std::min(vreco2->closest_vertex_distance_z, distance);
974  }
975  }
976  return recopv;
977 }
978 
979 // ------------ method called to produce the data ------------
980 void Primary4DVertexValidation::matchReco2Sim(std::vector<recoPrimaryVertex>& recopv,
981  std::vector<simPrimaryVertex>& simpv,
982  const edm::ValueMap<float>& sigmat0,
983  const edm::ValueMap<float>& MVA,
984  const edm::Handle<reco::BeamSpot>& BS) {
985  for (auto vv : simpv) {
986  vv.wnt.clear();
987  vv.wos.clear();
988  }
989  for (auto rv : recopv) {
990  rv.wnt.clear();
991  rv.wos.clear();
992  }
993 
994  for (unsigned int iv = 0; iv < recopv.size(); iv++) {
995  const reco::Vertex* vertex = recopv.at(iv).recVtx;
996 
997  for (unsigned int iev = 0; iev < simpv.size(); iev++) {
998  double wnt = 0;
999  double wos = 0;
1000  double evwnt = 0;
1001  double evwos = 0;
1002  double evnt = 0;
1003 
1004  for (auto iTrack = vertex->tracks_begin(); iTrack != vertex->tracks_end(); ++iTrack) {
1005  double pt = (*iTrack)->pt();
1006 
1007  if (vertex->trackWeight(*iTrack) < trackweightTh_)
1008  continue;
1009  if (MVA[(*iTrack)] < mvaTh_)
1010  continue;
1011 
1012  auto tp_info = getMatchedTP(*iTrack, simpv.at(iev).sim_vertex);
1013  if (tp_info != nullptr) {
1014  double dz2_beam = pow((*BS).BeamWidthX() * cos((*iTrack)->phi()) / tan((*iTrack)->theta()), 2) +
1015  pow((*BS).BeamWidthY() * sin((*iTrack)->phi()) / tan((*iTrack)->theta()), 2);
1016  double dz2 = pow((*iTrack)->dzError(), 2) + dz2_beam +
1017  pow(0.0020, 2); // added 20 um, some tracks have crazy small resolutions
1018  wos = vertex->trackWeight(*iTrack) / dz2;
1019  wnt = vertex->trackWeight(*iTrack) * std::min(pt, 1.0);
1020 
1021  if (sigmat0[(*iTrack)] > 0) {
1022  double sigmaZ = (*BS).sigmaZ();
1023  double sigmaT = sigmaZ / c_; // c in cm/ns
1024  wos = wos / erf(sigmat0[(*iTrack)] / sigmaT);
1025  }
1026  simpv.at(iev).addTrack(iv, wos, wnt);
1027  recopv.at(iv).addTrack(iev, wos, wnt);
1028  evwos += wos;
1029  evwnt += wnt;
1030  evnt++;
1031  }
1032  } //RecoTracks loop
1033 
1034  // require 2 tracks for a wos-match
1035  if ((evwos > 0) && (evwos > recopv.at(iv).maxwos) && (evnt > 1)) {
1036  recopv.at(iv).wosmatch = iev;
1037  recopv.at(iv).maxwos = evwos;
1038  recopv.at(iv).maxwosnt = evnt;
1039 
1040  simpv.at(iev).wos_dominated_recv.push_back(iv);
1041  simpv.at(iev).nwosmatch++;
1042  }
1043 
1044  // weighted track counting match, require at least one track
1045  if ((evnt > 0) && (evwnt > recopv.at(iv).maxwnt)) {
1046  recopv.at(iv).wntmatch = iev;
1047  recopv.at(iv).maxwnt = evwnt;
1048  }
1049  } //TrackingVertex loop
1050 
1051  } //RecoPrimaryVertex
1052 
1053  //after filling infos, goes for the sim-reco match
1054  for (auto& vrec : recopv) {
1055  vrec.sim = NOT_MATCHED;
1056  vrec.matchQuality = 0;
1057  }
1058  unsigned int iev = 0;
1059  for (auto& vv : simpv) {
1060  if (debug_) {
1061  edm::LogPrint("Primary4DVertexValidation") << "iev: " << iev;
1062  edm::LogPrint("Primary4DVertexValidation") << "wos_dominated_recv.size: " << vv.wos_dominated_recv.size();
1063  }
1064  for (unsigned int i = 0; i < vv.wos_dominated_recv.size(); i++) {
1065  auto recov = vv.wos_dominated_recv.at(i);
1066  if (debug_) {
1067  edm::LogPrint("Primary4DVertexValidation")
1068  << "index of reco vertex: " << recov << " that has a wos: " << vv.wos.at(recov) << " at position " << i;
1069  }
1070  }
1071  vv.rec = NOT_MATCHED;
1072  vv.matchQuality = 0;
1073  iev++;
1074  }
1075  //this tries a one-to-one match, taking simPV with highest wos if there are > 1 simPV candidates
1076  for (unsigned int rank = 1; rank < maxRank_; rank++) {
1077  for (unsigned int iev = 0; iev < simpv.size(); iev++) { //loop on SimPV
1078  if (simpv.at(iev).rec != NOT_MATCHED)
1079  continue;
1080  if (simpv.at(iev).nwosmatch == 0)
1081  continue;
1082  if (simpv.at(iev).nwosmatch > rank)
1083  continue;
1084  unsigned int iv = NOT_MATCHED;
1085  for (unsigned int k = 0; k < simpv.at(iev).wos_dominated_recv.size(); k++) {
1086  unsigned int rec = simpv.at(iev).wos_dominated_recv.at(k);
1087  auto vrec = recopv.at(rec);
1088  if (vrec.sim != NOT_MATCHED)
1089  continue; // already matched
1090  if (std::abs(simpv.at(iev).z - vrec.z) > zWosMatchMax_)
1091  continue; // insanely far away
1092  if ((iv == NOT_MATCHED) || simpv.at(iev).wos.at(rec) > simpv.at(iev).wos.at(iv)) {
1093  iv = rec;
1094  }
1095  }
1096  if (iv !=
1097  NOT_MATCHED) { //if the rec vertex has already been associated is possible that iv remains NOT_MATCHED at this point
1098  recopv.at(iv).sim = iev;
1099  simpv.at(iev).rec = iv;
1100  recopv.at(iv).matchQuality = rank;
1101  simpv.at(iev).matchQuality = rank;
1102  }
1103  }
1104  }
1105  //give vertices a chance that have a lot of overlap, but are still recognizably
1106  //caused by a specific simvertex (without being classified as dominating)
1107  //like a small peak sitting on the flank of a larger nearby peak
1108  unsigned int ntry = 0;
1109  while (ntry++ < maxTry_) {
1110  unsigned nmatch = 0;
1111  for (unsigned int iev = 0; iev < simpv.size(); iev++) {
1112  if ((simpv.at(iev).rec != NOT_MATCHED) || (simpv.at(iev).wos.empty()))
1113  continue;
1114  // find a rec vertex for the NOT_MATCHED sim vertex
1115  unsigned int rec = NOT_MATCHED;
1116  for (auto rv : simpv.at(iev).wos) {
1117  if ((rec == NOT_MATCHED) || (rv.second > simpv.at(iev).wos.at(rec))) {
1118  rec = rv.first;
1119  }
1120  }
1121 
1122  if (rec == NOT_MATCHED) { //try with wnt match
1123  for (auto rv : simpv.at(iev).wnt) {
1124  if ((rec == NOT_MATCHED) || (rv.second > simpv.at(iev).wnt.at(rec))) {
1125  rec = rv.first;
1126  }
1127  }
1128  }
1129 
1130  if (rec == NOT_MATCHED)
1131  continue;
1132  if (recopv.at(rec).sim != NOT_MATCHED)
1133  continue; // already gone
1134 
1135  // check if the recvertex can be matched
1136  unsigned int rec2sim = NOT_MATCHED;
1137  for (auto sv : recopv.at(rec).wos) {
1138  if (simpv.at(sv.first).rec != NOT_MATCHED)
1139  continue; // already used
1140  if ((rec2sim == NOT_MATCHED) || (sv.second > recopv.at(rec).wos.at(rec2sim))) {
1141  rec2sim = sv.first;
1142  }
1143  }
1144  if (iev == rec2sim) {
1145  // do the match and assign lowest quality (i.e. max rank)
1146  recopv.at(rec).sim = iev;
1147  recopv.at(rec).matchQuality = maxRank_;
1148  simpv.at(iev).rec = rec;
1149  simpv.at(iev).matchQuality = maxRank_;
1150  nmatch++;
1151  }
1152  } //sim loop
1153  if (nmatch == 0) {
1154  break;
1155  }
1156  } // ntry
1157 }
1158 
1159 void Primary4DVertexValidation::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
1160  using edm::Handle;
1161  using edm::View;
1162  using std::cout;
1163  using std::endl;
1164  using std::vector;
1165  using namespace reco;
1166 
1167  std::vector<float> pileUpInfo_z;
1168 
1169  // get the pileup information
1170  edm::Handle<std::vector<PileupSummaryInfo>> puinfoH;
1171  if (iEvent.getByToken(vecPileupSummaryInfoToken_, puinfoH)) {
1172  for (auto const& pu_info : *puinfoH.product()) {
1173  if (pu_info.getBunchCrossing() == 0) {
1174  pileUpInfo_z = pu_info.getPU_zpositions();
1175  break;
1176  }
1177  }
1178  }
1179 
1180  edm::Handle<TrackingParticleCollection> TPCollectionH;
1181  iEvent.getByToken(trackingParticleCollectionToken_, TPCollectionH);
1182  if (!TPCollectionH.isValid())
1183  edm::LogWarning("Primary4DVertexValidation") << "TPCollectionH is not valid";
1184 
1185  edm::Handle<TrackingVertexCollection> TVCollectionH;
1186  iEvent.getByToken(trackingVertexCollectionToken_, TVCollectionH);
1187  if (!TVCollectionH.isValid())
1188  edm::LogWarning("Primary4DVertexValidation") << "TVCollectionH is not valid";
1189 
1190  edm::Handle<reco::SimToRecoCollection> simToRecoH;
1191  iEvent.getByToken(simToRecoAssociationToken_, simToRecoH);
1192  if (simToRecoH.isValid())
1193  s2r_ = simToRecoH.product();
1194  else
1195  edm::LogWarning("Primary4DVertexValidation") << "simToRecoH is not valid";
1196 
1197  edm::Handle<reco::RecoToSimCollection> recoToSimH;
1198  iEvent.getByToken(recoToSimAssociationToken_, recoToSimH);
1199  if (recoToSimH.isValid())
1200  r2s_ = recoToSimH.product();
1201  else
1202  edm::LogWarning("Primary4DVertexValidation") << "recoToSimH is not valid";
1203 
1204  edm::Handle<reco::BeamSpot> BeamSpotH;
1205  iEvent.getByToken(RecBeamSpotToken_, BeamSpotH);
1206  if (!BeamSpotH.isValid())
1207  edm::LogWarning("Primary4DVertexValidation") << "BeamSpotH is not valid";
1208 
1209  std::vector<simPrimaryVertex> simpv; // a list of simulated primary MC vertices
1210  simpv = getSimPVs(TVCollectionH);
1211  //this bool check if first vertex in that with highest pT
1212  bool signal_is_highest_pt =
1213  std::max_element(simpv.begin(), simpv.end(), [](const simPrimaryVertex& lhs, const simPrimaryVertex& rhs) {
1214  return lhs.ptsq < rhs.ptsq;
1215  }) == simpv.begin();
1216 
1217  std::vector<recoPrimaryVertex> recopv; // a list of reconstructed primary MC vertices
1218  edm::Handle<edm::View<reco::Vertex>> recVtxs;
1219  iEvent.getByToken(Rec4DVerToken_, recVtxs);
1220  if (!recVtxs.isValid())
1221  edm::LogWarning("Primary4DVertexValidation") << "recVtxs is not valid";
1222  recopv = getRecoPVs(recVtxs);
1223 
1224  const auto& trackAssoc = iEvent.get(trackAssocToken_);
1225  const auto& pathLength = iEvent.get(pathLengthToken_);
1226  const auto& momentum = iEvent.get(momentumToken_);
1227  const auto& time = iEvent.get(timeToken_);
1228  const auto& t0Pid = iEvent.get(t0PidToken_);
1229  const auto& t0Safe = iEvent.get(t0SafePidToken_);
1230  const auto& sigmat0Safe = iEvent.get(sigmat0SafePidToken_);
1231  const auto& mtdQualMVA = iEvent.get(trackMVAQualToken_);
1232  const auto& tMtd = iEvent.get(tmtdToken_);
1233  const auto& tofPi = iEvent.get(tofPiToken_);
1234  const auto& tofK = iEvent.get(tofKToken_);
1235  const auto& tofP = iEvent.get(tofPToken_);
1236  const auto& probPi = iEvent.get(probPiToken_);
1237  const auto& probK = iEvent.get(probKToken_);
1238  const auto& probP = iEvent.get(probPToken_);
1239 
1240  //I have simPV and recoPV collections
1241  matchReco2Sim(recopv, simpv, sigmat0Safe, mtdQualMVA, BeamSpotH);
1242 
1243  //Loop on tracks
1244  for (unsigned int iv = 0; iv < recopv.size(); iv++) {
1245  if (recopv.at(iv).ndof > selNdof_) {
1246  const reco::Vertex* vertex = recopv.at(iv).recVtx;
1247 
1248  for (unsigned int iev = 0; iev < simpv.size(); iev++) {
1249  auto vsim = simpv.at(iev).sim_vertex;
1250 
1251  bool selectedVtxMatching = recopv.at(iv).sim == iev && simpv.at(iev).rec == iv &&
1252  simpv.at(iev).eventId.bunchCrossing() == 0 && simpv.at(iev).eventId.event() == 0 &&
1253  recopv.at(iv).OriginalIndex == 0;
1254  if (selectedVtxMatching && !recopv.at(iv).is_signal()) {
1255  edm::LogWarning("Primary4DVertexValidation")
1256  << "Reco vtx leading match inconsistent: BX/ID " << simpv.at(iev).eventId.bunchCrossing() << " "
1257  << simpv.at(iev).eventId.event();
1258  }
1259  double vzsim = simpv.at(iev).z;
1260  double vtsim = simpv.at(iev).t * simUnit_;
1261 
1262  for (auto iTrack = vertex->tracks_begin(); iTrack != vertex->tracks_end(); ++iTrack) {
1263  if (trackAssoc[*iTrack] == -1) {
1264  LogTrace("mtdTracks") << "Extended track not associated";
1265  continue;
1266  }
1267 
1268  if (vertex->trackWeight(*iTrack) < trackweightTh_)
1269  continue;
1270 
1271  bool noCrack = std::abs((*iTrack)->eta()) < trackMaxBtlEta_ || std::abs((*iTrack)->eta()) > trackMinEtlEta_;
1272 
1273  bool selectRecoTrk = mvaRecSel(**iTrack, *vertex, t0Safe[*iTrack], sigmat0Safe[*iTrack]);
1274  if (selectedVtxMatching && selectRecoTrk) {
1275  if (noCrack) {
1276  meMVATrackEffPtTot_->Fill((*iTrack)->pt());
1277  }
1278  meMVATrackEffEtaTot_->Fill(std::abs((*iTrack)->eta()));
1279  }
1280 
1281  auto tp_info = getMatchedTP(*iTrack, vsim);
1282  if (tp_info != nullptr) {
1283  double mass = (*tp_info)->mass();
1284  double tsim = (*tp_info)->parentVertex()->position().t() * simUnit_;
1285  double tEst = timeFromTrueMass(mass, pathLength[*iTrack], momentum[*iTrack], time[*iTrack]);
1286 
1287  double xsim = (*tp_info)->parentVertex()->position().x();
1288  double ysim = (*tp_info)->parentVertex()->position().y();
1289  double zsim = (*tp_info)->parentVertex()->position().z();
1290  double xPCA = (*iTrack)->vx();
1291  double yPCA = (*iTrack)->vy();
1292  double zPCA = (*iTrack)->vz();
1293 
1294  double dZ = zPCA - zsim;
1295  double d3D = std::sqrt((xPCA - xsim) * (xPCA - xsim) + (yPCA - ysim) * (yPCA - ysim) + dZ * dZ);
1296  // orient d3D according to the projection of RECO - SIM onto simulated momentum
1297  if ((xPCA - xsim) * ((*tp_info)->px()) + (yPCA - ysim) * ((*tp_info)->py()) + dZ * ((*tp_info)->pz()) <
1298  0.) {
1299  d3D = -d3D;
1300  }
1301 
1302  bool selectTP = mvaTPSel(**tp_info);
1303 
1304  if (selectedVtxMatching && selectRecoTrk && selectTP) {
1305  meMVATrackZposResTot_->Fill((*iTrack)->vz() - vzsim);
1306  if (noCrack) {
1307  meMVATrackMatchedEffPtTot_->Fill((*iTrack)->pt());
1308  }
1309  meMVATrackMatchedEffEtaTot_->Fill(std::abs((*iTrack)->eta()));
1310  }
1311 
1312  if (sigmat0Safe[*iTrack] == -1)
1313  continue;
1314 
1315  if (selectedVtxMatching && selectRecoTrk && selectTP) {
1316  meMVATrackResTot_->Fill(t0Safe[*iTrack] - vtsim);
1317  meMVATrackPullTot_->Fill((t0Safe[*iTrack] - vtsim) / sigmat0Safe[*iTrack]);
1318  if (noCrack) {
1319  meMVATrackMatchedEffPtMtd_->Fill((*iTrack)->pt());
1320  }
1321  meMVATrackMatchedEffEtaMtd_->Fill(std::abs((*iTrack)->eta()));
1322 
1323  unsigned int noPIDtype = 0;
1324  if (probPi[*iTrack] == -1) {
1325  noPIDtype = 1;
1326  } else if (std::isnan(probPi[*iTrack])) {
1327  noPIDtype = 2;
1328  } else if (probPi[*iTrack] == 1 && probK[*iTrack] == 0 && probP[*iTrack] == 0) {
1329  noPIDtype = 3;
1330  }
1331  bool noPID = noPIDtype > 0;
1332  bool isPi = !noPID && 1. - probPi[*iTrack] < minProbHeavy_;
1333  bool isK = !noPID && !isPi && probK[*iTrack] > probP[*iTrack];
1334  bool isP = !noPID && !isPi && !isK;
1335 
1336  if ((isPi && std::abs(tMtd[*iTrack] - tofPi[*iTrack] - t0Pid[*iTrack]) > tol_) ||
1337  (isK && std::abs(tMtd[*iTrack] - tofK[*iTrack] - t0Pid[*iTrack]) > tol_) ||
1338  (isP && std::abs(tMtd[*iTrack] - tofP[*iTrack] - t0Pid[*iTrack]) > tol_)) {
1339  edm::LogWarning("Primary4DVertexValidation")
1340  << "No match between mass hyp. and time: " << std::abs((*tp_info)->pdgId()) << " mass hyp pi/k/p "
1341  << isPi << " " << isK << " " << isP << " t0/t0safe " << t0Pid[*iTrack] << " " << t0Safe[*iTrack]
1342  << " tMtd - tof pi/K/p " << tMtd[*iTrack] - tofPi[*iTrack] << " " << tMtd[*iTrack] - tofK[*iTrack]
1343  << " " << tMtd[*iTrack] - tofP[*iTrack] << " Prob pi/K/p " << probPi[*iTrack] << " "
1344  << probK[*iTrack] << " " << probP[*iTrack];
1345  }
1346 
1347  if (std::abs((*iTrack)->eta()) < trackMaxBtlEta_) {
1348  meBarrelPIDp_->Fill((*iTrack)->p());
1349  meBarrelNoPIDtype_->Fill(noPIDtype + 0.5);
1350  if (std::abs((*tp_info)->pdgId()) == 211) {
1351  if (noPID) {
1352  meBarrelTruePiNoPID_->Fill((*iTrack)->p());
1353  } else if (isPi) {
1354  meBarrelTruePiAsPi_->Fill((*iTrack)->p());
1355  } else if (isK) {
1356  meBarrelTruePiAsK_->Fill((*iTrack)->p());
1357  } else if (isP) {
1358  meBarrelTruePiAsP_->Fill((*iTrack)->p());
1359  } else {
1360  edm::LogWarning("Primary4DVertexValidation")
1361  << "No PID class: " << std::abs((*tp_info)->pdgId()) << " t0/t0safe " << t0Pid[*iTrack] << " "
1362  << t0Safe[*iTrack] << " Prob pi/K/p " << probPi[*iTrack] << " " << probK[*iTrack] << " "
1363  << probP[*iTrack];
1364  }
1365  } else if (std::abs((*tp_info)->pdgId()) == 321) {
1366  if (noPID) {
1367  meBarrelTrueKNoPID_->Fill((*iTrack)->p());
1368  } else if (isPi) {
1369  meBarrelTrueKAsPi_->Fill((*iTrack)->p());
1370  } else if (isK) {
1371  meBarrelTrueKAsK_->Fill((*iTrack)->p());
1372  } else if (isP) {
1373  meBarrelTrueKAsP_->Fill((*iTrack)->p());
1374  } else {
1375  edm::LogWarning("Primary4DVertexValidation")
1376  << "No PID class: " << std::abs((*tp_info)->pdgId()) << " t0/t0safe " << t0Pid[*iTrack] << " "
1377  << t0Safe[*iTrack] << " Prob pi/K/p " << probPi[*iTrack] << " " << probK[*iTrack] << " "
1378  << probP[*iTrack];
1379  }
1380  } else if (std::abs((*tp_info)->pdgId()) == 2212) {
1381  if (noPID) {
1382  meBarrelTruePNoPID_->Fill((*iTrack)->p());
1383  } else if (isPi) {
1384  meBarrelTruePAsPi_->Fill((*iTrack)->p());
1385  } else if (isK) {
1386  meBarrelTruePAsK_->Fill((*iTrack)->p());
1387  } else if (isP) {
1388  meBarrelTruePAsP_->Fill((*iTrack)->p());
1389  } else {
1390  edm::LogWarning("Primary4DVertexValidation")
1391  << "No PID class: " << std::abs((*tp_info)->pdgId()) << " t0/t0safe " << t0Pid[*iTrack] << " "
1392  << t0Safe[*iTrack] << " Prob pi/K/p " << probPi[*iTrack] << " " << probK[*iTrack] << " "
1393  << probP[*iTrack];
1394  }
1395  }
1396  } else if (std::abs((*iTrack)->eta()) > trackMinEtlEta_ && std::abs((*iTrack)->eta()) < trackMaxEtlEta_) {
1397  meEndcapPIDp_->Fill((*iTrack)->p());
1398  meEndcapNoPIDtype_->Fill(noPIDtype + 0.5);
1399  if (std::abs((*tp_info)->pdgId()) == 211) {
1400  if (noPID) {
1401  meEndcapTruePiNoPID_->Fill((*iTrack)->p());
1402  } else if (isPi) {
1403  meEndcapTruePiAsPi_->Fill((*iTrack)->p());
1404  } else if (isK) {
1405  meEndcapTruePiAsK_->Fill((*iTrack)->p());
1406  } else if (isP) {
1407  meEndcapTruePiAsP_->Fill((*iTrack)->p());
1408  } else {
1409  edm::LogWarning("Primary4DVertexValidation")
1410  << "No PID class: " << std::abs((*tp_info)->pdgId()) << " t0/t0safe " << t0Pid[*iTrack] << " "
1411  << t0Safe[*iTrack] << " Prob pi/K/p " << probPi[*iTrack] << " " << probK[*iTrack] << " "
1412  << probP[*iTrack];
1413  }
1414  } else if (std::abs((*tp_info)->pdgId()) == 321) {
1415  if (noPID) {
1416  meEndcapTrueKNoPID_->Fill((*iTrack)->p());
1417  } else if (isPi) {
1418  meEndcapTrueKAsPi_->Fill((*iTrack)->p());
1419  } else if (isK) {
1420  meEndcapTrueKAsK_->Fill((*iTrack)->p());
1421  } else if (isP) {
1422  meEndcapTrueKAsP_->Fill((*iTrack)->p());
1423  } else {
1424  edm::LogWarning("Primary4DVertexValidation")
1425  << "No PID class: " << std::abs((*tp_info)->pdgId()) << " t0/t0safe " << t0Pid[*iTrack] << " "
1426  << t0Safe[*iTrack] << " Prob pi/K/p " << probPi[*iTrack] << " " << probK[*iTrack] << " "
1427  << probP[*iTrack];
1428  }
1429  } else if (std::abs((*tp_info)->pdgId()) == 2212) {
1430  if (noPID) {
1431  meEndcapTruePNoPID_->Fill((*iTrack)->p());
1432  } else if (isPi) {
1433  meEndcapTruePAsPi_->Fill((*iTrack)->p());
1434  } else if (isK) {
1435  meEndcapTruePAsK_->Fill((*iTrack)->p());
1436  } else if (isP) {
1437  meEndcapTruePAsP_->Fill((*iTrack)->p());
1438  } else {
1439  edm::LogWarning("Primary4DVertexValidation")
1440  << "No PID class: " << std::abs((*tp_info)->pdgId()) << " t0/t0safe " << t0Pid[*iTrack] << " "
1441  << t0Safe[*iTrack] << " Prob pi/K/p " << probPi[*iTrack] << " " << probK[*iTrack] << " "
1442  << probP[*iTrack];
1443  }
1444  }
1445  }
1446  }
1447  meTrackResTot_->Fill(t0Safe[*iTrack] - tsim);
1448  meTrackPullTot_->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
1449  meTrackZposResTot_->Fill(dZ);
1450  if ((*iTrack)->p() <= 2) {
1451  meTrackResLowPTot_->Fill(t0Safe[*iTrack] - tsim);
1452  meTrackPullLowPTot_->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
1453  } else {
1454  meTrackResHighPTot_->Fill(t0Safe[*iTrack] - tsim);
1455  meTrackPullHighPTot_->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
1456  }
1457 
1458  if (mtdQualMVA[(*iTrack)] < mvaL_) {
1459  meTrackZposRes_[0]->Fill(dZ);
1460  meTrack3DposRes_[0]->Fill(d3D);
1461  meTrackRes_[0]->Fill(t0Safe[*iTrack] - tsim);
1462  meTrackPull_[0]->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
1463 
1464  if (optionalPlots_) {
1465  meTrackResMass_[0]->Fill(t0Safe[*iTrack] - tEst);
1466  meTrackResMassTrue_[0]->Fill(tEst - tsim);
1467  }
1468 
1469  if ((*iTrack)->p() <= 2) {
1470  meTrackResLowP_[0]->Fill(t0Safe[*iTrack] - tsim);
1471  meTrackPullLowP_[0]->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
1472  } else if ((*iTrack)->p() > 2) {
1473  meTrackResHighP_[0]->Fill(t0Safe[*iTrack] - tsim);
1474  meTrackPullHighP_[0]->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
1475  }
1476 
1477  if (optionalPlots_) {
1478  if (std::abs((*tp_info)->pdgId()) == 2212) {
1479  meTrackResMassProtons_[0]->Fill(t0Safe[*iTrack] - tEst);
1480  meTrackResMassTrueProtons_[0]->Fill(tEst - tsim);
1481  } else if (std::abs((*tp_info)->pdgId()) == 211) {
1482  meTrackResMassPions_[0]->Fill(t0Safe[*iTrack] - tEst);
1483  meTrackResMassTruePions_[0]->Fill(tEst - tsim);
1484  }
1485  }
1486 
1487  } else if (mtdQualMVA[(*iTrack)] > mvaL_ && mtdQualMVA[(*iTrack)] < mvaH_) {
1488  meTrackZposRes_[1]->Fill(dZ);
1489  meTrack3DposRes_[1]->Fill(d3D);
1490  meTrackRes_[1]->Fill(t0Safe[*iTrack] - tsim);
1491  meTrackPull_[1]->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
1492 
1493  if (optionalPlots_) {
1494  meTrackResMass_[1]->Fill(t0Safe[*iTrack] - tEst);
1495  meTrackResMassTrue_[1]->Fill(tEst - tsim);
1496  }
1497 
1498  if ((*iTrack)->p() <= 2) {
1499  meTrackResLowP_[1]->Fill(t0Safe[*iTrack] - tsim);
1500  meTrackPullLowP_[1]->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
1501  } else if ((*iTrack)->p() > 2) {
1502  meTrackResHighP_[1]->Fill(t0Safe[*iTrack] - tsim);
1503  meTrackPullHighP_[1]->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
1504  }
1505 
1506  if (optionalPlots_) {
1507  if (std::abs((*tp_info)->pdgId()) == 2212) {
1508  meTrackResMassProtons_[1]->Fill(t0Safe[*iTrack] - tEst);
1509  meTrackResMassTrueProtons_[1]->Fill(tEst - tsim);
1510  } else if (std::abs((*tp_info)->pdgId()) == 211) {
1511  meTrackResMassPions_[1]->Fill(t0Safe[*iTrack] - tEst);
1512  meTrackResMassTruePions_[1]->Fill(tEst - tsim);
1513  }
1514  }
1515 
1516  } else if (mtdQualMVA[(*iTrack)] > mvaH_) {
1517  meTrackZposRes_[2]->Fill(dZ);
1518  meTrack3DposRes_[2]->Fill(d3D);
1519  meTrackRes_[2]->Fill(t0Safe[*iTrack] - tsim);
1520  meTrackPull_[2]->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
1521 
1522  if (optionalPlots_) {
1523  meTrackResMass_[2]->Fill(t0Safe[*iTrack] - tEst);
1524  meTrackResMassTrue_[2]->Fill(tEst - tsim);
1525  }
1526 
1527  if ((*iTrack)->p() <= 2) {
1528  meTrackResLowP_[2]->Fill(t0Safe[*iTrack] - tsim);
1529  meTrackPullLowP_[2]->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
1530  } else if ((*iTrack)->p() > 2) {
1531  meTrackResHighP_[2]->Fill(t0Safe[*iTrack] - tsim);
1532  meTrackPullHighP_[2]->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
1533  }
1534 
1535  if (optionalPlots_) {
1536  if (std::abs((*tp_info)->pdgId()) == 2212) {
1537  meTrackResMassProtons_[2]->Fill(t0Safe[*iTrack] - tEst);
1538  meTrackResMassTrueProtons_[2]->Fill(tEst - tsim);
1539  } else if (std::abs((*tp_info)->pdgId()) == 211) {
1540  meTrackResMassPions_[2]->Fill(t0Safe[*iTrack] - tEst);
1541  meTrackResMassTruePions_[2]->Fill(tEst - tsim);
1542  }
1543  }
1544  }
1545  } //if tp_info != nullptr
1546  }
1547  }
1548  } // ndof
1549  }
1550 
1551  int real = 0;
1552  int fake = 0;
1553  int other_fake = 0;
1554  int split = 0;
1555 
1556  auto puLineDensity = [&](double z) {
1557  // gaussian parameterization of line density vs z, z in cm, parameters in mm
1558  double argl = (z * 10. - lineDensityPar_[1]) / lineDensityPar_[2];
1559  return lineDensityPar_[0] * exp(-0.5 * argl * argl);
1560  };
1561 
1562  meRecVerNumber_->Fill(recopv.size());
1563  for (unsigned int ir = 0; ir < recopv.size(); ir++) {
1564  if (recopv.at(ir).ndof > selNdof_) {
1565  meRecoVtxVsLineDensity_->Fill(puLineDensity(recopv.at(ir).z));
1566  meRecPVZ_->Fill(recopv.at(ir).z, 1. / puLineDensity(recopv.at(ir).z));
1567  if (recopv.at(ir).recVtx->tError() > 0.) {
1568  meRecPVT_->Fill(recopv.at(ir).recVtx->t());
1569  }
1570  if (debug_) {
1571  edm::LogPrint("Primary4DVertexValidation") << "************* IR: " << ir;
1572  edm::LogPrint("Primary4DVertexValidation")
1573  << "z: " << recopv.at(ir).z << " corresponding to line density: " << puLineDensity(recopv.at(ir).z);
1574  edm::LogPrint("Primary4DVertexValidation") << "is_real: " << recopv.at(ir).is_real();
1575  edm::LogPrint("Primary4DVertexValidation") << "is_fake: " << recopv.at(ir).is_fake();
1576  edm::LogPrint("Primary4DVertexValidation") << "is_signal: " << recopv.at(ir).is_signal();
1577  edm::LogPrint("Primary4DVertexValidation") << "split_from: " << recopv.at(ir).split_from();
1578  edm::LogPrint("Primary4DVertexValidation") << "other fake: " << recopv.at(ir).other_fake();
1579  }
1580  if (recopv.at(ir).is_real())
1581  real++;
1582  if (recopv.at(ir).is_fake())
1583  fake++;
1584  if (recopv.at(ir).other_fake())
1585  other_fake++;
1586  if (recopv.at(ir).split_from() != -1) {
1587  split++;
1588  }
1589  } // ndof
1590  }
1591 
1592  if (debug_) {
1593  edm::LogPrint("Primary4DVertexValidation") << "is_real: " << real;
1594  edm::LogPrint("Primary4DVertexValidation") << "is_fake: " << fake;
1595  edm::LogPrint("Primary4DVertexValidation") << "split_from: " << split;
1596  edm::LogPrint("Primary4DVertexValidation") << "other fake: " << other_fake;
1597  }
1598  mePUvsRealV_->Fill(simpv.size(), real);
1599  mePUvsOtherFakeV_->Fill(simpv.size(), other_fake);
1600  mePUvsSplitV_->Fill(simpv.size(), split);
1601 
1602  //fill vertices histograms here in a new loop
1603  for (unsigned int is = 0; is < simpv.size(); is++) {
1604  // protect against particle guns with very displaced vertices
1605  if (std::isinf(1. / puLineDensity(simpv.at(is).z))) {
1606  continue;
1607  }
1608  meSimPVZ_->Fill(simpv.at(is).z, 1. / puLineDensity(simpv.at(is).z));
1609  if (is == 0 && optionalPlots_) {
1610  meSimPosInSimOrigCollection_->Fill(simpv.at(is).OriginalIndex);
1611  }
1612 
1613  if (simpv.at(is).rec == NOT_MATCHED) {
1614  if (debug_) {
1615  edm::LogPrint("Primary4DVertexValidation") << "sim vertex: " << is << " is not matched with any reco";
1616  }
1617  continue;
1618  }
1619 
1620  for (unsigned int ir = 0; ir < recopv.size(); ir++) {
1621  if (recopv.at(ir).ndof > selNdof_) {
1622  if (recopv.at(ir).sim == is && simpv.at(is).rec == ir) {
1623  meTimeRes_->Fill(recopv.at(ir).recVtx->t() - simpv.at(is).t * simUnit_);
1624  meTimePull_->Fill((recopv.at(ir).recVtx->t() - simpv.at(is).t * simUnit_) / recopv.at(ir).recVtx->tError());
1625  meMatchQual_->Fill(recopv.at(ir).matchQuality - 0.5);
1626  if (ir == 0) { //signal vertex plots
1627  meTimeSignalRes_->Fill(recopv.at(ir).recVtx->t() - simpv.at(is).t * simUnit_);
1628  meTimeSignalPull_->Fill((recopv.at(ir).recVtx->t() - simpv.at(is).t * simUnit_) /
1629  recopv.at(ir).recVtx->tError());
1630  if (optionalPlots_) {
1631  meRecoPosInSimCollection_->Fill(recopv.at(ir).sim);
1632  meRecoPosInRecoOrigCollection_->Fill(recopv.at(ir).OriginalIndex);
1633  }
1634  }
1635  if (simpv.at(is).eventId.bunchCrossing() == 0 && simpv.at(is).eventId.event() == 0) {
1636  if (!recopv.at(ir).is_signal()) {
1637  edm::LogWarning("Primary4DVertexValidation")
1638  << "Reco vtx leading match inconsistent: BX/ID " << simpv.at(is).eventId.bunchCrossing() << " "
1639  << simpv.at(is).eventId.event();
1640  }
1641  meRecoPVPosSignal_->Fill(
1642  recopv.at(ir).OriginalIndex); // position in reco vtx correction associated to sim signal
1643  if (!signal_is_highest_pt) {
1644  meRecoPVPosSignalNotHighestPt_->Fill(
1645  recopv.at(ir).OriginalIndex); // position in reco vtx correction associated to sim signal
1646  }
1647  }
1648 
1649  if (debug_) {
1650  edm::LogPrint("Primary4DVertexValidation") << "*** Matching RECO: " << ir << "with SIM: " << is << " ***";
1651  edm::LogPrint("Primary4DVertexValidation") << "Match Quality is " << recopv.at(ir).matchQuality;
1652  edm::LogPrint("Primary4DVertexValidation") << "****";
1653  }
1654  }
1655  } // ndof
1656  }
1657  }
1658 
1659  //dz histos
1660  for (unsigned int iv = 0; iv < recVtxs->size() - 1; iv++) {
1661  if (recVtxs->at(iv).ndof() > selNdof_) {
1662  double mindistance_realreal = 1e10;
1663 
1664  for (unsigned int jv = iv; jv < recVtxs->size(); jv++) {
1665  if ((!(jv == iv)) && select(recVtxs->at(jv))) {
1666  double dz = recVtxs->at(iv).z() - recVtxs->at(jv).z();
1667  double dtsigma = std::sqrt(recVtxs->at(iv).covariance(3, 3) + recVtxs->at(jv).covariance(3, 3));
1668  double dt = (std::abs(dz) <= deltaZcut_ && dtsigma > 0.)
1669  ? (recVtxs->at(iv).t() - recVtxs->at(jv).t()) / dtsigma
1670  : -9999.;
1671  if (recopv.at(iv).is_real() && recopv.at(jv).is_real()) {
1672  meDeltaZrealreal_->Fill(std::abs(dz));
1673  if (dt != -9999.) {
1674  meDeltaTrealreal_->Fill(std::abs(dt));
1675  }
1676  if (std::abs(dz) < std::abs(mindistance_realreal)) {
1677  mindistance_realreal = dz;
1678  }
1679  } else if (recopv.at(iv).is_fake() && recopv.at(jv).is_fake()) {
1680  meDeltaZfakefake_->Fill(std::abs(dz));
1681  if (dt != -9999.) {
1682  meDeltaTfakefake_->Fill(std::abs(dt));
1683  }
1684  }
1685  }
1686  }
1687 
1688  double mindistance_fakereal = 1e10;
1689  double mindistance_realfake = 1e10;
1690  for (unsigned int jv = 0; jv < recVtxs->size(); jv++) {
1691  if ((!(jv == iv)) && select(recVtxs->at(jv))) {
1692  double dz = recVtxs->at(iv).z() - recVtxs->at(jv).z();
1693  double dtsigma = std::sqrt(recVtxs->at(iv).covariance(3, 3) + recVtxs->at(jv).covariance(3, 3));
1694  double dt = (std::abs(dz) <= deltaZcut_ && dtsigma > 0.)
1695  ? (recVtxs->at(iv).t() - recVtxs->at(jv).t()) / dtsigma
1696  : -9999.;
1697  if (recopv.at(iv).is_fake() && recopv.at(jv).is_real()) {
1698  meDeltaZfakereal_->Fill(std::abs(dz));
1699  if (dt != -9999.) {
1700  meDeltaTfakereal_->Fill(std::abs(dt));
1701  }
1702  if (std::abs(dz) < std::abs(mindistance_fakereal)) {
1703  mindistance_fakereal = dz;
1704  }
1705  }
1706 
1707  if (recopv.at(iv).is_real() && recopv.at(jv).is_fake() && (std::abs(dz) < std::abs(mindistance_realfake))) {
1708  mindistance_realfake = dz;
1709  }
1710  }
1711  }
1712  } //ndof
1713  }
1714 
1715 } // end of analyze
1716 
1717 void Primary4DVertexValidation::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
1718  edm::ParameterSetDescription desc;
1719 
1720  desc.add<std::string>("folder", "MTD/Vertices");
1721  desc.add<edm::InputTag>("TPtoRecoTrackAssoc", edm::InputTag("trackingParticleRecoTrackAsssociation"));
1722  desc.add<edm::InputTag>("mtdTracks", edm::InputTag("trackExtenderWithMTD"));
1723  desc.add<edm::InputTag>("SimTag", edm::InputTag("mix", "MergedTrackTruth"));
1724  desc.add<edm::InputTag>("offlineBS", edm::InputTag("offlineBeamSpot"));
1725  desc.add<edm::InputTag>("offline4DPV", edm::InputTag("offlinePrimaryVertices4D"));
1726  desc.add<edm::InputTag>("trackAssocSrc", edm::InputTag("trackExtenderWithMTD:generalTrackassoc"))
1727  ->setComment("Association between General and MTD Extended tracks");
1728  desc.add<edm::InputTag>("pathLengthSrc", edm::InputTag("trackExtenderWithMTD:generalTrackPathLength"));
1729  desc.add<edm::InputTag>("momentumSrc", edm::InputTag("trackExtenderWithMTD:generalTrackp"));
1730  desc.add<edm::InputTag>("tmtd", edm::InputTag("trackExtenderWithMTD:generalTracktmtd"));
1731  desc.add<edm::InputTag>("timeSrc", edm::InputTag("trackExtenderWithMTD:generalTracktmtd"));
1732  desc.add<edm::InputTag>("sigmaSrc", edm::InputTag("trackExtenderWithMTD:generalTracksigmatmtd"));
1733  desc.add<edm::InputTag>("t0PID", edm::InputTag("tofPID:t0"));
1734  desc.add<edm::InputTag>("t0SafePID", edm::InputTag("tofPID:t0safe"));
1735  desc.add<edm::InputTag>("sigmat0SafePID", edm::InputTag("tofPID:sigmat0safe"));
1736  desc.add<edm::InputTag>("trackMVAQual", edm::InputTag("mtdTrackQualityMVA:mtdQualMVA"));
1737  desc.add<edm::InputTag>("tofPi", edm::InputTag("trackExtenderWithMTD:generalTrackTofPi"));
1738  desc.add<edm::InputTag>("tofK", edm::InputTag("trackExtenderWithMTD:generalTrackTofK"));
1739  desc.add<edm::InputTag>("tofP", edm::InputTag("trackExtenderWithMTD:generalTrackTofP"));
1740  desc.add<edm::InputTag>("probPi", edm::InputTag("tofPID:probPi"));
1741  desc.add<edm::InputTag>("probK", edm::InputTag("tofPID:probK"));
1742  desc.add<edm::InputTag>("probP", edm::InputTag("tofPID:probP"));
1743  desc.add<bool>("useOnlyChargedTracks", true);
1744  desc.addUntracked<bool>("debug", false);
1745  desc.addUntracked<bool>("optionalPlots", false);
1746  desc.add<double>("trackweightTh", 0.5);
1747  desc.add<double>("mvaTh", 0.01);
1748  desc.add<double>("minProbHeavy", 0.75);
1749 
1750  //lineDensity parameters have been obtained by fitting the distribution of the z position of the vertices,
1751  //using a 200k single mu ptGun sample (gaussian fit)
1752  std::vector<double> lDP;
1753  lDP.push_back(1.87);
1754  lDP.push_back(0.);
1755  lDP.push_back(42.5);
1756  desc.add<std::vector<double>>("lineDensityPar", lDP);
1757  descriptions.add("vertices4DValid", desc);
1758 }
1759 
1760 const bool Primary4DVertexValidation::mvaTPSel(const TrackingParticle& tp) {
1761  bool match = false;
1762  if (tp.status() != 1) {
1763  return match;
1764  }
1765  match = tp.charge() != 0 && tp.pt() > pTcut_ && std::abs(tp.eta()) < etacutGEN_;
1766  return match;
1767 }
1768 
1769 const bool Primary4DVertexValidation::mvaRecSel(const reco::TrackBase& trk,
1770  const reco::Vertex& vtx,
1771  const double& t0,
1772  const double& st0) {
1773  bool match = false;
1774  match = trk.pt() > pTcut_ && std::abs(trk.eta()) < etacutREC_ && std::abs(trk.vz() - vtx.z()) <= deltaZcut_;
1775  if (st0 > 0.) {
1776  match = match && std::abs(t0 - vtx.t()) < 3. * st0;
1777  }
1778  return match;
1779 }
1780 
1781 DEFINE_FWK_MODULE(Primary4DVertexValidation);
FrontierCondition_GT_autoExpress_cfi.t0
t0
Definition: FrontierCondition_GT_autoExpress_cfi.py:149
Primary4DVertexValidation::trackMaxBtlEta_
static constexpr double trackMaxBtlEta_
Definition: Primary4DVertexValidation.cc:226
dt
float dt
Definition: AMPTWrapper.h:136
Primary4DVertexValidation::Rec4DVerToken_
edm::EDGetTokenT< edm::View< reco::Vertex > > Rec4DVerToken_
Definition: Primary4DVertexValidation.cc:253
CommonMethods.isnan
def isnan(num)
Definition: CommonMethods.py:97
ProducerED_cfi.InputTag
InputTag
Definition: ProducerED_cfi.py:5
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
Definition: ParameterSet.h:307
Primary4DVertexValidation::trackingVertexCollectionToken_
edm::EDGetTokenT< TrackingVertexCollection > trackingVertexCollectionToken_
Definition: Primary4DVertexValidation.cc:249
bphysicsOniaDQM_cfi.vertex
vertex
Definition: bphysicsOniaDQM_cfi.py:7
Primary4DVertexValidation::mvaTPSel
const bool mvaTPSel(const TrackingParticle &)
Definition: Primary4DVertexValidation.cc:1760
Primary4DVertexValidation::trackMaxEtlEta_
static constexpr double trackMaxEtlEta_
Definition: Primary4DVertexValidation.cc:228
mps_fire.i
i
Definition: mps_fire.py:429
Primary4DVertexValidation::etacutREC_
static constexpr double etacutREC_
Definition: Primary4DVertexValidation.cc:223
folder_
std::string folder_
Definition: DQMEDAnalyzer.cc:60
muonTagProbeFilters_cff.matched
matched
Definition: muonTagProbeFilters_cff.py:62
Primary4DVertexValidation::trackingParticleCollectionToken_
edm::EDGetTokenT< TrackingParticleCollection > trackingParticleCollectionToken_
Definition: Primary4DVertexValidation.cc:248
dqm::implementation::NavigatorBase::setCurrentFolder
virtual void setCurrentFolder(std::string const &fullpath)
Definition: DQMStore.cc:36
Primary4DVertexValidation::sigmat0SafePidToken_
edm::EDGetTokenT< edm::ValueMap< float > > sigmat0SafePidToken_
Definition: Primary4DVertexValidation.cc:262
Primary4DVertexValidation::getSimPVs
std::vector< Primary4DVertexValidation::simPrimaryVertex > getSimPVs(const edm::Handle< TrackingVertexCollection > &)
Definition: Primary4DVertexValidation.cc:815
Primary4DVertexValidation::RecBeamSpotToken_
edm::EDGetTokenT< reco::BeamSpot > RecBeamSpotToken_
Definition: Primary4DVertexValidation.cc:252
Primary4DVertexValidation::select
bool select(const reco::Vertex &, int level=0)
Definition: Primary4DVertexValidation.cc:793
Primary4DVertexValidation::matchReco2Sim
void matchReco2Sim(std::vector< recoPrimaryVertex > &, std::vector< simPrimaryVertex > &, const edm::ValueMap< float > &, const edm::ValueMap< float > &, const edm::Handle< reco::BeamSpot > &)
Definition: Primary4DVertexValidation.cc:980
Primary4DVertexValidation::probKToken_
edm::EDGetTokenT< edm::ValueMap< float > > probKToken_
Definition: Primary4DVertexValidation.cc:270
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
reco::VertexBaseRef
edm::RefToBase< reco::Vertex > VertexBaseRef
persistent reference to a Vertex, using views
Definition: VertexFwd.h:21
edm::Handle::product
T const * product() const
Definition: Handle.h:70
L1BJetProducer_cff.vtx
vtx
Definition: L1BJetProducer_cff.py:10
Primary4DVertexValidation::matchRecoTrack2SimSignal
bool matchRecoTrack2SimSignal(const reco::TrackBaseRef &)
Definition: Primary4DVertexValidation.cc:744
caHitNtupletGeneratorKernels::iev
TkSoAView< TrackerTraits > HitToTuple< TrackerTraits > const *__restrict__ int32_t int32_t int iev
Definition: CAHitNtupletGeneratorKernelsImpl.h:858
Primary4DVertexValidation::etacutGEN_
static constexpr double etacutGEN_
Definition: Primary4DVertexValidation.cc:222
l1ctLayer1_cff.dZ
dZ
Definition: l1ctLayer1_cff.py:69
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
pfDeepBoostedJetPreprocessParams_cfi.sv
sv
Definition: pfDeepBoostedJetPreprocessParams_cfi.py:352
ALPAKA_ACCELERATOR_NAMESPACE::brokenline::constexpr
if constexpr(n > 3)
Definition: BrokenLine.h:164
cms::cuda::assert
assert(be >=bs)
Primary4DVertexValidation::analyze
void analyze(const edm::Event &, const edm::EventSetup &) override
Definition: Primary4DVertexValidation.cc:1159
Primary4DVertexValidation::timeToken_
edm::EDGetTokenT< edm::ValueMap< float > > timeToken_
Definition: Primary4DVertexValidation.cc:258
Primary4DVertexValidation::getRecoPVs
std::vector< Primary4DVertexValidation::recoPrimaryVertex > getRecoPVs(const edm::Handle< edm::View< reco::Vertex >> &)
Definition: Primary4DVertexValidation.cc:916
hcalRecHitTable_cff.time
time
Definition: hcalRecHitTable_cff.py:14
Primary4DVertexValidation::trackMinEtlEta_
static constexpr double trackMinEtlEta_
Definition: Primary4DVertexValidation.cc:227
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
edm::AssociationMap::find
const_iterator find(const key_type &k) const
find element with specified reference key
Definition: AssociationMap.h:175
LogTrace
#define LogTrace(id)
Definition: MessageLogger.h:242
DiDispStaMuonMonitor_cfi.pt
pt
Definition: DiDispStaMuonMonitor_cfi.py:39
Primary4DVertexValidation::meTrackResMassProtons_
MonitorElement * meTrackResMassProtons_[3]
Definition: Primary4DVertexValidation.cc:328
Primary4DVertexValidation::meTrackResMassTrueProtons_
MonitorElement * meTrackResMassTrueProtons_[3]
Definition: Primary4DVertexValidation.cc:329
edm::AssociationMap::end
const_iterator end() const
last iterator over the map (read only)
Definition: AssociationMap.h:173
dqm::impl::MonitorElement::Fill
void Fill(long long x)
Definition: MonitorElement.h:292
math::XYZTLorentzVector
XYZTLorentzVectorD XYZTLorentzVector
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:29
reco::TrackBase::pt
double pt() const
track transverse momentum
Definition: TrackBase.h:637
Primary4DVertexValidation::timeFromTrueMass
double timeFromTrueMass(double, double, double, double)
Definition: Primary4DVertexValidation.cc:781
iEvent
int iEvent
Definition: GenABIO.cc:224
Primary4DVertexValidation::LorentzVector
math::XYZTLorentzVector LorentzVector
Definition: Primary4DVertexValidation.cc:55
ALPAKA_ACCELERATOR_NAMESPACE::vertexFinder::iv
int32_t *__restrict__ iv
Definition: clusterTracksByDensity.h:58
reco::TrackBase::vz
double vz() const
z coordinate of the reference point on track
Definition: TrackBase.h:661
Primary4DVertexValidation::simPrimaryVertex::addTrack
void addTrack(unsigned int irecv, double twos, double twt)
Definition: Primary4DVertexValidation.cc:99
dqm::implementation::IBooker::bookProfile
MonitorElement * bookProfile(TString const &name, TString const &title, int nchX, double lowX, double highX, int, double lowY, double highY, char const *option="s", FUNC onbooking=NOOP())
Definition: DQMStore.h:408
V0Monitor_cff.v0
v0
Definition: V0Monitor_cff.py:7
Primary4DVertexValidation::fillDescriptions
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
Definition: Primary4DVertexValidation.cc:1717
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
Primary4DVertexValidation::Primary4DVertexValidation
Primary4DVertexValidation(const edm::ParameterSet &)
Definition: Primary4DVertexValidation.cc:369
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
Primary4DVertexValidation::probPiToken_
edm::EDGetTokenT< edm::ValueMap< float > > probPiToken_
Definition: Primary4DVertexValidation.cc:269
funct::tan
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
Primary4DVertexValidation::probPToken_
edm::EDGetTokenT< edm::ValueMap< float > > probPToken_
Definition: Primary4DVertexValidation.cc:271
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
Primary4DVertexValidation::momentumToken_
edm::EDGetTokenT< edm::ValueMap< float > > momentumToken_
Definition: Primary4DVertexValidation.cc:257
Primary4DVertexValidation::meTrackResMassTruePions_
MonitorElement * meTrackResMassTruePions_[3]
Definition: Primary4DVertexValidation.cc:331
Primary4DVertexValidation::tofKToken_
edm::EDGetTokenT< edm::ValueMap< float > > tofKToken_
Definition: Primary4DVertexValidation.cc:267
Primary4DVertexValidation::RecTrackToken_
edm::EDGetTokenT< reco::TrackCollection > RecTrackToken_
Definition: Primary4DVertexValidation.cc:244
HLT_2024v14_cff.distance
distance
Definition: HLT_2024v14_cff.py:5740
Primary4DVertexValidation::lineDensityPar_
const std::vector< double > lineDensityPar_
Definition: Primary4DVertexValidation.cc:240
edm::LogPrint
Log< level::Warning, true > LogPrint
Definition: MessageLogger.h:136
Primary4DVertexValidation::r2s_
const reco::RecoToSimCollection * r2s_
Definition: Primary4DVertexValidation.cc:241
Primary4DVertexValidation::vecPileupSummaryInfoToken_
edm::EDGetTokenT< std::vector< PileupSummaryInfo > > vecPileupSummaryInfoToken_
Definition: Primary4DVertexValidation.cc:246
Primary4DVertexValidation::simPrimaryVertex::simPrimaryVertex
simPrimaryVertex(double x1, double y1, double z1, double t1)
Definition: Primary4DVertexValidation.cc:59
Primary4DVertexValidation::deltaZcut_
static constexpr double deltaZcut_
Definition: Primary4DVertexValidation.cc:225
EgHLTOffHistBins_cfi.mass
mass
Definition: EgHLTOffHistBins_cfi.py:34
Primary4DVertexValidation::mvaRecSel
const bool mvaRecSel(const reco::TrackBase &, const reco::Vertex &, const double &, const double &)
Definition: Primary4DVertexValidation.cc:1769
Primary4DVertexValidation::t0PidToken_
edm::EDGetTokenT< edm::ValueMap< float > > t0PidToken_
Definition: Primary4DVertexValidation.cc:260
TrackingVertexRef
edm::Ref< TrackingVertexCollection > TrackingVertexRef
Definition: TrackingVertexContainer.h:9
Primary4DVertexValidation::NOT_MATCHED
static constexpr unsigned int NOT_MATCHED
Definition: Primary4DVertexValidation.cc:213
reco::TrackBase::eta
double eta() const
pseudorapidity of momentum vector
Definition: TrackBase.h:652
Primary4DVertexValidation::simToRecoAssociationToken_
edm::EDGetTokenT< reco::SimToRecoCollection > simToRecoAssociationToken_
Definition: Primary4DVertexValidation.cc:250
Primary4DVertexValidation::trackAssocToken_
edm::EDGetTokenT< edm::ValueMap< int > > trackAssocToken_
Definition: Primary4DVertexValidation.cc:255
angle_units::operators::convertMmToCm
constexpr NumType convertMmToCm(NumType millimeters)
Definition: angle_units.h:44
Primary4DVertexValidation::getMatchedTP
const edm::Ref< std::vector< TrackingParticle > > * getMatchedTP(const reco::TrackBaseRef &, const TrackingVertexRef &)
Definition: Primary4DVertexValidation.cc:761
edm::ConfigurationDescriptions::add
void add(std::string const &label, ParameterSetDescription const &psetDescription)
Definition: ConfigurationDescriptions.cc:57
nano_mu_local_reco_cff.bool
bool
Definition: nano_mu_local_reco_cff.py:14
testProducerWithPsetDescEmpty_cfi.y1
y1
Definition: testProducerWithPsetDescEmpty_cfi.py:29
Primary4DVertexValidation::recoPrimaryVertex::addTrack
void addTrack(unsigned int iev, double twos, double wt)
Definition: Primary4DVertexValidation.cc:169
edm::HandleBase::isValid
bool isValid() const
Definition: HandleBase.h:70
reco
fixed size matrix
Definition: AlignmentAlgorithmBase.h:46
Primary4DVertexValidation::s2r_
const reco::SimToRecoCollection * s2r_
Definition: Primary4DVertexValidation.cc:242
Primary4DVertexValidation::tofPiToken_
edm::EDGetTokenT< edm::ValueMap< float > > tofPiToken_
Definition: Primary4DVertexValidation.cc:266
TrackingParticle
Monte Carlo truth information used for tracking validation.
Definition: TrackingParticle.h:29
Primary4DVertexValidation::bookHistograms
void bookHistograms(DQMStore::IBooker &i, edm::Run const &, edm::EventSetup const &) override
Definition: Primary4DVertexValidation.cc:411
Primary4DVertexValidation::recoToSimAssociationToken_
edm::EDGetTokenT< reco::RecoToSimCollection > recoToSimAssociationToken_
Definition: Primary4DVertexValidation.cc:251
match
std::pair< typename Association::data_type::first_type, double > match(Reference key, Association association, bool bestMatchByMaxValue)
Generic matching function.
Definition: Utils.h:10
gather_cfg.cout
cout
Definition: gather_cfg.py:144
JetChargeProducer_cfi.exp
exp
Definition: JetChargeProducer_cfi.py:6
edm::LogWarning
Log< level::Warning, false > LogWarning
Definition: MessageLogger.h:128
dqm::implementation::IBooker::book1D
MonitorElement * book1D(TString const &name, TString const &title, int const nchX, double const lowX, double const highX, FUNC onbooking=NOOP())
Definition: DQMStore.h:98
Primary4DVertexValidation::t0SafePidToken_
edm::EDGetTokenT< edm::ValueMap< float > > t0SafePidToken_
Definition: Primary4DVertexValidation.cc:261
Primary4DVertexValidation::tmtdToken_
edm::EDGetTokenT< edm::ValueMap< float > > tmtdToken_
Definition: Primary4DVertexValidation.cc:265
testProducerWithPsetDescEmpty_cfi.x1
x1
Definition: testProducerWithPsetDescEmpty_cfi.py:33
Primary4DVertexValidation::zWosMatchMax_
static constexpr double zWosMatchMax_
Definition: Primary4DVertexValidation.cc:221
Primary4DVertexValidation::tofPToken_
edm::EDGetTokenT< edm::ValueMap< float > > tofPToken_
Definition: Primary4DVertexValidation.cc:268
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
edm::Run
Definition: Run.h:45
Primary4DVertexValidation::trackMVAQualToken_
edm::EDGetTokenT< edm::ValueMap< float > > trackMVAQualToken_
Definition: Primary4DVertexValidation.cc:263
Primary4DVertexValidation::pathLengthToken_
edm::EDGetTokenT< edm::ValueMap< float > > pathLengthToken_
Definition: Primary4DVertexValidation.cc:256
Generated for CMSSW Reference Manual by doxygen 1.8.14