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