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