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