CMS 3D CMS Logo
Primary4DVertexValidation.cc
Go to the documentation of this file.
1 #include <numeric>
2 #include <string>
3 #include <vector>
4 #include <map>
5 #include <algorithm>
6 
7 #include "FWCore/Framework/interface/Frameworkfwd.h"
8 #include "FWCore/Framework/interface/ESHandle.h"
9 #include "FWCore/Utilities/interface/isFinite.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 #include "FWCore/Utilities/interface/isFinite.h"
18 
19 #include "DataFormats/Math/interface/LorentzVector.h"
20 #include "DataFormats/Math/interface/Point3D.h"
21 
22 // reco track and vertex
23 #include "DataFormats/TrackReco/interface/Track.h"
24 #include "DataFormats/TrackReco/interface/TrackFwd.h"
25 #include "DataFormats/VertexReco/interface/Vertex.h"
26 #include "DataFormats/VertexReco/interface/VertexFwd.h"
27 #include "RecoVertex/VertexPrimitives/interface/TransientVertex.h"
28 
29 // TrackingParticle
30 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
31 #include "SimDataFormats/Associations/interface/TrackToTrackingParticleAssociator.h"
32 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h"
33 #include "SimDataFormats/TrackingAnalysis/interface/TrackingVertexContainer.h"
34 
35 // Fastjet
36 #include <fastjet/internal/base.hh>
37 #include "fastjet/PseudoJet.hh"
38 #include "fastjet/JetDefinition.hh"
39 #include "fastjet/ClusterSequence.hh"
40 #include "fastjet/Selector.hh"
41 
42 // HepPDTRecord
43 #include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"
44 
45 // pile-up
46 #include "SimDataFormats/PileupSummaryInfo/interface/PileupSummaryInfo.h"
47 
48 // associator
49 #include "SimTracker/VertexAssociation/interface/calculateVertexSharedTracks.h"
50 
51 // vertexing
52 #include "RecoVertex/PrimaryVertexProducer/interface/TrackFilterForPVFinding.h"
53 
54 // simulated vertex
55 #include "SimDataFormats/Associations/interface/VertexToTrackingVertexAssociator.h"
56 
57 // DQM
58 #include "DQMServices/Core/interface/DQMEDAnalyzer.h"
59 #include "DQMServices/Core/interface/DQMStore.h"
60 
61 #include "DataFormats/Math/interface/deltaR.h"
62 #include "DataFormats/Math/interface/GeantUnits.h"
63 #include "CLHEP/Units/PhysicalConstants.h"
64 
65 // class declaration
66 class Primary4DVertexValidation : public DQMEDAnalyzer {
67  typedef math::XYZTLorentzVector LorentzVector;
68 
69  // auxiliary class holding simulated vertices
70  struct simPrimaryVertex {
71  simPrimaryVertex(double x1, double y1, double z1, double t1)
72  : x(x1),
73  y(y1),
74  z(z1),
75  t(t1),
76  pt(0),
77  ptsq(0),
78  closest_vertex_distance_z(-1.),
79  nGenTrk(0),
80  num_matched_reco_tracks(0),
81  average_match_quality(0.0) {
82  ptot.setPx(0);
83  ptot.setPy(0);
84  ptot.setPz(0);
85  ptot.setE(0);
86  p4 = LorentzVector(0, 0, 0, 0);
87  r = sqrt(x * x + y * y);
88  };
89  double x, y, z, r, t;
90  HepMC::FourVector ptot;
91  LorentzVector p4;
92  double pt;
93  double ptsq;
94  double closest_vertex_distance_z;
95  int nGenTrk;
96  int num_matched_reco_tracks;
97  float average_match_quality;
98  EncodedEventId eventId;
99  TrackingVertexRef sim_vertex;
100  int OriginalIndex = -1;
101 
102  unsigned int nwosmatch = 0; // number of rec vertices dominated by this sim evt (by wos)
103  unsigned int nwntmatch = 0; // number of rec vertices dominated by this sim evt (by wnt)
104  std::vector<unsigned int> wos_dominated_recv; // list of dominated rec vertices (by wos, size==nwosmatch)
105 
106  std::map<unsigned int, double> wnt; // weighted number of tracks in recvtx (by index)
107  std::map<unsigned int, double> wos; // weight over sigma**2 in recvtx (by index)
108  double sumwos = 0; // sum of wos in any recvtx
109  double sumwnt = 0; // sum of weighted tracks
110  unsigned int rec = NOT_MATCHED; // best match (NO_MATCH if not matched)
111  unsigned int matchQuality = 0; // quality flag
112 
113  void addTrack(unsigned int irecv, double twos, double twt) {
114  sumwnt += twt;
115  if (wnt.find(irecv) == wnt.end()) {
116  wnt[irecv] = twt;
117  } else {
118  wnt[irecv] += twt;
119  }
120 
121  sumwos += twos;
122  if (wos.find(irecv) == wos.end()) {
123  wos[irecv] = twos;
124  } else {
125  wos[irecv] += twos;
126  }
127  }
128  };
129 
130  // auxiliary class holding reconstructed vertices
131  struct recoPrimaryVertex {
132  recoPrimaryVertex(double x1, double y1, double z1)
133  : x(x1),
134  y(y1),
135  z(z1),
136  pt(0),
137  ptsq(0),
138  closest_vertex_distance_z(-1.),
139  nRecoTrk(0),
140  num_matched_sim_tracks(0),
141  ndof(0.),
142  recVtx(nullptr) {
143  r = sqrt(x * x + y * y);
144  };
145  double x, y, z, r;
146  double pt;
147  double ptsq;
148  double closest_vertex_distance_z;
149  int nRecoTrk;
150  int num_matched_sim_tracks;
151  double ndof;
152  const reco::Vertex* recVtx;
153  reco::VertexBaseRef recVtxRef;
154  int OriginalIndex = -1;
155 
156  std::map<unsigned int, double> wos; // sim event -> wos
157  std::map<unsigned int, double> wnt; // sim event -> weighted number of truth matched tracks
158  unsigned int wosmatch = NOT_MATCHED; // index of the sim event providing the largest contribution to wos
159  unsigned int wntmatch = NOT_MATCHED; // index of the sim event providing the highest number of tracks
160  double sumwos = 0; // total sum of wos of all truth matched tracks
161  double sumwnt = 0; // total weighted number of truth matchted tracks
162  double maxwos = 0; // largest wos sum from one sim event (wosmatch)
163  double maxwnt = 0; // largest weighted number of tracks from one sim event (ntmatch)
164  unsigned int maxwosnt = 0; // number of tracks from the sim event with highest wos
165  unsigned int sim = NOT_MATCHED; // best match (NO_MATCH if not matched)
166  unsigned int matchQuality = 0; // quality flag
167 
168  bool is_real() { return (matchQuality > 0) && (matchQuality < 99); }
169 
170  bool is_fake() { return (matchQuality <= 0) || (matchQuality >= 99); }
171 
172  bool is_signal() { return (sim == 0); }
173 
174  int split_from() {
175  if (is_real())
176  return -1;
177  if ((maxwos > 0) && (maxwos > 0.3 * sumwos))
178  return wosmatch;
179  return -1;
180  }
181  bool other_fake() { return (is_fake() && (split_from() < 0)); }
182 
183  void addTrack(unsigned int iev, double twos, double wt) {
184  sumwnt += wt;
185  if (wnt.find(iev) == wnt.end()) {
186  wnt[iev] = wt;
187  } else {
188  wnt[iev] += wt;
189  }
190 
191  sumwos += twos;
192  if (wos.find(iev) == wos.end()) {
193  wos[iev] = twos;
194  } else {
195  wos[iev] += twos;
196  }
197  }
198  };
199 
200 public:
201  explicit Primary4DVertexValidation(const edm::ParameterSet&);
202  ~Primary4DVertexValidation() override;
203 
204  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
205  void analyze(const edm::Event&, const edm::EventSetup&) override;
206  void bookHistograms(DQMStore::IBooker& i, edm::Run const&, edm::EventSetup const&) override;
207 
208 private:
209  void matchReco2Sim(std::vector<recoPrimaryVertex>&,
210  std::vector<simPrimaryVertex>&,
211  const edm::ValueMap<float>&,
212  const edm::ValueMap<float>&,
213  const edm::Handle<reco::BeamSpot>&);
214  bool matchRecoTrack2SimSignal(const reco::TrackBaseRef&);
215  std::pair<const edm::Ref<std::vector<TrackingParticle>>*, int> getMatchedTP(const reco::TrackBaseRef&,
216  const TrackingVertexRef&);
217  double timeFromTrueMass(double, double, double, double);
218  bool select(const reco::Vertex&, int level = 0);
219  void observablesFromJets(const std::vector<reco::Track>&,
220  const std::vector<double>&,
221  const std::vector<int>&,
222  const std::string&,
223  unsigned int&,
224  double&,
225  double&,
226  double&,
227  double&);
228  void isParticle(const reco::TrackBaseRef&,
229  const edm::ValueMap<float>&,
230  const edm::ValueMap<float>&,
231  const edm::ValueMap<float>&,
232  const edm::ValueMap<float>&,
233  const edm::ValueMap<float>&,
234  unsigned int&,
235  bool&,
236  bool&,
237  bool&,
238  bool&);
239  void getWosWnt(const reco::Vertex&,
240  const reco::TrackBaseRef&,
241  const edm::ValueMap<float>&,
242  const edm::ValueMap<float>&,
243  const edm::Handle<reco::BeamSpot>&,
244  double&,
245  double&);
246  std::vector<Primary4DVertexValidation::simPrimaryVertex> getSimPVs(const edm::Handle<TrackingVertexCollection>&);
247  std::vector<Primary4DVertexValidation::recoPrimaryVertex> getRecoPVs(const edm::Handle<edm::View<reco::Vertex>>&);
248 
249  void printMatchedRecoTrackInfo(const reco::Vertex&,
250  const reco::TrackBaseRef&,
251  const TrackingParticleRef&,
252  const unsigned int&);
253  void printSimVtxRecoVtxInfo(const struct Primary4DVertexValidation::simPrimaryVertex&,
254  const struct Primary4DVertexValidation::recoPrimaryVertex&);
255  const bool trkTPSelLV(const TrackingParticle&);
256  const bool trkRecSel(const reco::TrackBase&);
257 
258  // ----------member data ---------------------------
259 
260  const std::string folder_;
261  static constexpr unsigned int NOT_MATCHED = 66666;
262  static constexpr double simUnit_ = 1e9; // sim time in s while reco time in ns
263  static constexpr double c_ = 2.99792458e1; // c in cm/ns
264  static constexpr double mvaL_ = 0.5; // MVA cuts for MVA categories
265  static constexpr double mvaH_ = 0.8;
266  static constexpr double selNdof_ = 4.;
267  static constexpr double maxRank_ = 8.;
268  static constexpr double maxTry_ = 10.;
269  static constexpr double zWosMatchMax_ = 1.;
270  static constexpr double etacutGEN_ = 4.; // |eta| < 4;
271  static constexpr double etacutREC_ = 3.; // |eta| < 3;
272  static constexpr double pTcut_ = 0.7; // PT > 0.7 GeV
273  static constexpr double deltaZcut_ = 0.1; // dz separation 1 mm
274  static constexpr double trackMaxBtlEta_ = 1.5;
275  static constexpr double trackMinEtlEta_ = 1.6;
276  static constexpr double trackMaxEtlEta_ = 3.;
277  static constexpr double tol_ = 1.e-4; // tolerance on reconstructed track time, [ns]
278  static constexpr double minThrSumWnt_ = 0.01; // min threshold for filling histograms with logarithmic scale
279  static constexpr double minThrSumWos_ = 0.1;
280  static constexpr double minThrSumPt_ = 0.01;
281  static constexpr double minThrSumPt2_ = 1.e-3;
282  static constexpr double minThrMetPt_ = 1.e-3;
283  static constexpr double minThrSumPz_ = 1.e-4;
284  static constexpr double rBTL_ = 110.0;
285  static constexpr double zETL_ = 290.0;
286 
287  static constexpr float c_cm_ns = geant_units::operators::convertMmToCm(CLHEP::c_light); // [mm/ns] -> [cm/ns]
288 
289  bool use_only_charged_tracks_;
290  bool optionalPlots_;
291  bool use3dNoTime_;
292 
293  const double minProbHeavy_;
294  const double trackweightTh_;
295  const double mvaTh_;
296  const std::vector<double> lineDensityPar_;
297  const reco::RecoToSimCollection* r2s_;
298  const reco::SimToRecoCollection* s2r_;
299 
300  edm::EDGetTokenT<reco::TrackCollection> RecTrackToken_;
301 
302  edm::EDGetTokenT<std::vector<PileupSummaryInfo>> vecPileupSummaryInfoToken_;
303 
304  edm::EDGetTokenT<TrackingParticleCollection> trackingParticleCollectionToken_;
305  edm::EDGetTokenT<TrackingVertexCollection> trackingVertexCollectionToken_;
306  edm::EDGetTokenT<reco::SimToRecoCollection> simToRecoAssociationToken_;
307  edm::EDGetTokenT<reco::RecoToSimCollection> recoToSimAssociationToken_;
308  edm::EDGetTokenT<reco::BeamSpot> RecBeamSpotToken_;
309  edm::EDGetTokenT<edm::View<reco::Vertex>> Rec4DVerToken_;
310 
311  edm::EDGetTokenT<edm::ValueMap<int>> trackAssocToken_;
312  edm::EDGetTokenT<edm::ValueMap<float>> pathLengthToken_;
313  edm::EDGetTokenT<edm::ValueMap<float>> momentumToken_;
314  edm::EDGetTokenT<edm::ValueMap<float>> timeToken_;
315 
316  edm::EDGetTokenT<edm::ValueMap<float>> t0PidToken_;
317  edm::EDGetTokenT<edm::ValueMap<float>> sigmat0PidToken_;
318  edm::EDGetTokenT<edm::ValueMap<float>> t0SafePidToken_;
319  edm::EDGetTokenT<edm::ValueMap<float>> sigmat0SafePidToken_;
320  edm::EDGetTokenT<edm::ValueMap<float>> trackMVAQualToken_;
321 
322  edm::EDGetTokenT<edm::ValueMap<float>> tmtdToken_;
323  edm::EDGetTokenT<edm::ValueMap<float>> tofPiToken_;
324  edm::EDGetTokenT<edm::ValueMap<float>> tofKToken_;
325  edm::EDGetTokenT<edm::ValueMap<float>> tofPToken_;
326  edm::EDGetTokenT<edm::ValueMap<float>> probPiToken_;
327  edm::EDGetTokenT<edm::ValueMap<float>> probKToken_;
328  edm::EDGetTokenT<edm::ValueMap<float>> probPToken_;
329  edm::ESGetToken<HepPDT::ParticleDataTable, edm::DefaultRecord> pdtToken_;
330 
331  // histogram declaration
332  MonitorElement* meTrackEffPtTot_;
333  MonitorElement* meTrackMatchedTPEffPtTot_;
334  MonitorElement* meTrackMatchedTPEffPtMtd_;
335  MonitorElement* meTrackEffEtaTot_;
336  MonitorElement* meTrackMatchedTPEffEtaTot_;
337  MonitorElement* meTrackMatchedTPEffEtaMtd_;
338  MonitorElement* meTrackMatchedTPResTot_;
339  MonitorElement* meTrackMatchedTPPullTot_;
340  MonitorElement* meTrackResTot_;
341  MonitorElement* meTrackPullTot_;
342  MonitorElement* meTrackRes_[3];
343  MonitorElement* meTrackPull_[3];
344  MonitorElement* meTrackResMass_[3];
345  MonitorElement* meTrackResMassTrue_[3];
346  MonitorElement* meTrackMatchedTPZposResTot_;
347  MonitorElement* meTrackZposResTot_;
348  MonitorElement* meTrackZposRes_[3];
349  MonitorElement* meTrack3DposRes_[3];
350  MonitorElement* meTimeRes_;
351  MonitorElement* meTimePull_;
352  MonitorElement* meTimeSignalRes_;
353  MonitorElement* meTimeSignalPull_;
354  MonitorElement* mePUvsRealV_;
355  MonitorElement* mePUvsFakeV_;
356  MonitorElement* mePUvsOtherFakeV_;
357  MonitorElement* mePUvsSplitV_;
358  MonitorElement* meMatchQual_;
359  MonitorElement* meDeltaZrealreal_;
360  MonitorElement* meDeltaZfakefake_;
361  MonitorElement* meDeltaZfakereal_;
362  MonitorElement* meDeltaTrealreal_;
363  MonitorElement* meDeltaTfakefake_;
364  MonitorElement* meDeltaTfakereal_;
365  MonitorElement* meRecoPosInSimCollection_;
366  MonitorElement* meRecoPosInRecoOrigCollection_;
367  MonitorElement* meSimPosInSimOrigCollection_;
368  MonitorElement* meRecoPVPosSignal_;
369  MonitorElement* meRecoPVPosSignalNotHighestPt_;
370  MonitorElement* meRecoVtxVsLineDensity_;
371  MonitorElement* meRecVerNumber_;
372  MonitorElement* meRecPVZ_;
373  MonitorElement* meRecPVT_;
374  MonitorElement* meSimPVZ_;
375 
376  MonitorElement* meVtxTrackMult_;
377  MonitorElement* meVtxTrackW_;
378  MonitorElement* meVtxTrackWnt_;
379  MonitorElement* meVtxTrackRecLVMult_;
380  MonitorElement* meVtxTrackRecLVW_;
381  MonitorElement* meVtxTrackRecLVWnt_;
382 
383  MonitorElement* mePUTrackMult_;
384  MonitorElement* mePUTrackRelMult_;
385  MonitorElement* meFakeTrackRelMult_;
386  MonitorElement* mePUTrackSumWnt_;
387  MonitorElement* mePUTrackRelSumWnt_;
388  MonitorElement* mePUTrackSumWos_;
389  MonitorElement* mePUTrackRelSumWos_;
390  MonitorElement* meSecTrackSumWos_;
391  MonitorElement* meSecTrackRelSumWos_;
392  MonitorElement* meFakeTrackRelSumWos_;
393  MonitorElement* mePUTrackWnt_;
394  MonitorElement* mePUTrackSumPt_;
395  MonitorElement* mePUTrackRelSumPt_;
396  MonitorElement* mePUTrackSumPt2_;
397  MonitorElement* mePUTrackRelSumPt2_;
398 
399  MonitorElement* mePUTrackRelMultvsMult_;
400  MonitorElement* meFakeTrackRelMultvsMult_;
401  MonitorElement* mePUTrackRelSumWntvsSumWnt_;
402  MonitorElement* mePUTrackRelSumWosvsSumWos_;
403  MonitorElement* meFakeTrackRelSumWosvsSumWos_;
404  MonitorElement* meSecTrackRelSumWosvsSumWos_;
405  MonitorElement* mePUTrackRelSumPtvsSumPt_;
406  MonitorElement* mePUTrackRelSumPt2vsSumPt2_;
407 
408  MonitorElement* mePUTrackRecLVMult_;
409  MonitorElement* mePUTrackRecLVRelMult_;
410  MonitorElement* meFakeTrackRecLVRelMult_;
411  MonitorElement* mePUTrackRecLVSumWnt_;
412  MonitorElement* mePUTrackRecLVRelSumWnt_;
413  MonitorElement* mePUTrackRecLVSumWos_;
414  MonitorElement* mePUTrackRecLVRelSumWos_;
415  MonitorElement* meSecTrackRecLVSumWos_;
416  MonitorElement* meSecTrackRecLVRelSumWos_;
417  MonitorElement* meFakeTrackRecLVRelSumWos_;
418  MonitorElement* mePUTrackRecLVWnt_;
419  MonitorElement* mePUTrackRecLVSumPt_;
420  MonitorElement* mePUTrackRecLVRelSumPt_;
421  MonitorElement* mePUTrackRecLVSumPt2_;
422  MonitorElement* mePUTrackRecLVRelSumPt2_;
423 
424  MonitorElement* mePUTrackRecLVRelMultvsMult_;
425  MonitorElement* meFakeTrackRecLVRelMultvsMult_;
426  MonitorElement* mePUTrackRecLVRelSumWntvsSumWnt_;
427  MonitorElement* mePUTrackRecLVRelSumWosvsSumWos_;
428  MonitorElement* meSecTrackRecLVRelSumWosvsSumWos_;
429  MonitorElement* meFakeTrackRecLVRelSumWosvsSumWos_;
430  MonitorElement* mePUTrackRecLVRelSumPtvsSumPt_;
431  MonitorElement* mePUTrackRecLVRelSumPt2vsSumPt2_;
432 
433  MonitorElement* meJetsPUMult_;
434  MonitorElement* meJetsPUHt_;
435  MonitorElement* meJetsPUSumPt2_;
436  MonitorElement* meJetsPUMetPt_;
437  MonitorElement* meJetsPUSumPz_;
438  MonitorElement* meJetsPURelMult_;
439  MonitorElement* meJetsPURelHt_;
440  MonitorElement* meJetsPURelSumPt2_;
441  MonitorElement* meJetsFakeRelSumPt2_;
442  MonitorElement* meJetsPURelMetPt_;
443  MonitorElement* meJetsPURelSumPz_;
444 
445  MonitorElement* meJetsPURelMultvsMult_;
446  MonitorElement* meJetsPURelHtvsHt_;
447  MonitorElement* meJetsPURelSumPt2vsSumPt2_;
448  MonitorElement* meJetsFakeRelSumPt2vsSumPt2_;
449  MonitorElement* meJetsPURelMetPtvsMetPt_;
450  MonitorElement* meJetsPURelSumPzvsSumPz_;
451 
452  MonitorElement* meJetsRecLVPUMult_;
453  MonitorElement* meJetsRecLVPUHt_;
454  MonitorElement* meJetsRecLVPUSumPt2_;
455  MonitorElement* meJetsRecLVPUMetPt_;
456  MonitorElement* meJetsRecLVPUSumPz_;
457  MonitorElement* meJetsRecLVPURelMult_;
458  MonitorElement* meJetsRecLVPURelHt_;
459  MonitorElement* meJetsRecLVPURelSumPt2_;
460  MonitorElement* meJetsRecLVFakeRelSumPt2_;
461  MonitorElement* meJetsRecLVPURelMetPt_;
462  MonitorElement* meJetsRecLVPURelSumPz_;
463 
464  MonitorElement* meJetsRecLVPURelMultvsMult_;
465  MonitorElement* meJetsRecLVPURelHtvsHt_;
466  MonitorElement* meJetsRecLVPURelSumPt2vsSumPt2_;
467  MonitorElement* meJetsRecLVFakeRelSumPt2vsSumPt2_;
468  MonitorElement* meJetsRecLVPURelMetPtvsMetPt_;
469  MonitorElement* meJetsRecLVPURelSumPzvsSumPz_;
470 
471  // some tests
472  MonitorElement* meTrackResLowPTot_;
473  MonitorElement* meTrackResHighPTot_;
474  MonitorElement* meTrackPullLowPTot_;
475  MonitorElement* meTrackPullHighPTot_;
476 
477  MonitorElement* meTrackResLowP_[3];
478  MonitorElement* meTrackResHighP_[3];
479  MonitorElement* meTrackPullLowP_[3];
480  MonitorElement* meTrackPullHighP_[3];
481 
482  MonitorElement* meTrackResMassProtons_[3];
483  MonitorElement* meTrackResMassTrueProtons_[3];
484  MonitorElement* meTrackResMassPions_[3];
485  MonitorElement* meTrackResMassTruePions_[3];
486 
487  MonitorElement* meBarrelPIDp_;
488  MonitorElement* meEndcapPIDp_;
489 
490  MonitorElement* meBarrelNoPIDtype_;
491  MonitorElement* meEndcapNoPIDtype_;
492 
493  MonitorElement* meBarrelTruePiNoPID_;
494  MonitorElement* meBarrelTrueKNoPID_;
495  MonitorElement* meBarrelTruePNoPID_;
496  MonitorElement* meEndcapTruePiNoPID_;
497  MonitorElement* meEndcapTrueKNoPID_;
498  MonitorElement* meEndcapTruePNoPID_;
499 
500  MonitorElement* meBarrelTruePiAsPi_;
501  MonitorElement* meBarrelTruePiAsK_;
502  MonitorElement* meBarrelTruePiAsP_;
503  MonitorElement* meEndcapTruePiAsPi_;
504  MonitorElement* meEndcapTruePiAsK_;
505  MonitorElement* meEndcapTruePiAsP_;
506 
507  MonitorElement* meBarrelTrueKAsPi_;
508  MonitorElement* meBarrelTrueKAsK_;
509  MonitorElement* meBarrelTrueKAsP_;
510  MonitorElement* meEndcapTrueKAsPi_;
511  MonitorElement* meEndcapTrueKAsK_;
512  MonitorElement* meEndcapTrueKAsP_;
513 
514  MonitorElement* meBarrelTruePAsPi_;
515  MonitorElement* meBarrelTruePAsK_;
516  MonitorElement* meBarrelTruePAsP_;
517  MonitorElement* meEndcapTruePAsPi_;
518  MonitorElement* meEndcapTruePAsK_;
519  MonitorElement* meEndcapTruePAsP_;
520 };
521 
522 // constructors and destructor
523 Primary4DVertexValidation::Primary4DVertexValidation(const edm::ParameterSet& iConfig)
524  : folder_(iConfig.getParameter<std::string>("folder")),
525  use_only_charged_tracks_(iConfig.getParameter<bool>("useOnlyChargedTracks")),
526  optionalPlots_(iConfig.getUntrackedParameter<bool>("optionalPlots")),
527  use3dNoTime_(iConfig.getParameter<bool>("use3dNoTime")),
528  minProbHeavy_(iConfig.getParameter<double>("minProbHeavy")),
529  trackweightTh_(iConfig.getParameter<double>("trackweightTh")),
530  mvaTh_(iConfig.getParameter<double>("mvaTh")),
531  lineDensityPar_(iConfig.getParameter<std::vector<double>>("lineDensityPar")),
532  pdtToken_(esConsumes<HepPDT::ParticleDataTable, edm::DefaultRecord>()) {
533  vecPileupSummaryInfoToken_ = consumes<std::vector<PileupSummaryInfo>>(edm::InputTag(std::string("addPileupInfo")));
534  trackingParticleCollectionToken_ =
535  consumes<TrackingParticleCollection>(iConfig.getParameter<edm::InputTag>("SimTag"));
536  trackingVertexCollectionToken_ = consumes<TrackingVertexCollection>(iConfig.getParameter<edm::InputTag>("SimTag"));
537  simToRecoAssociationToken_ =
538  consumes<reco::SimToRecoCollection>(iConfig.getParameter<edm::InputTag>("TPtoRecoTrackAssoc"));
539  recoToSimAssociationToken_ =
540  consumes<reco::RecoToSimCollection>(iConfig.getParameter<edm::InputTag>("TPtoRecoTrackAssoc"));
541  RecTrackToken_ = consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("mtdTracks"));
542  RecBeamSpotToken_ = consumes<reco::BeamSpot>(iConfig.getParameter<edm::InputTag>("offlineBS"));
543  Rec4DVerToken_ = consumes<edm::View<reco::Vertex>>(iConfig.getParameter<edm::InputTag>("offline4DPV"));
544  trackAssocToken_ = consumes<edm::ValueMap<int>>(iConfig.getParameter<edm::InputTag>("trackAssocSrc"));
545  pathLengthToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("pathLengthSrc"));
546  momentumToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("momentumSrc"));
547  timeToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("timeSrc"));
548  t0PidToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("t0PID"));
549  sigmat0PidToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("sigmat0PID"));
550  t0SafePidToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("t0SafePID"));
551  sigmat0SafePidToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("sigmat0SafePID"));
552  trackMVAQualToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("trackMVAQual"));
553  tmtdToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("tmtd"));
554  tofPiToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("tofPi"));
555  tofKToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("tofK"));
556  tofPToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("tofP"));
557  probPiToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("probPi"));
558  probKToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("probK"));
559  probPToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("probP"));
560 }
561 
562 Primary4DVertexValidation::~Primary4DVertexValidation() {}
563 
564 //
565 // member functions
566 //
567 void Primary4DVertexValidation::bookHistograms(DQMStore::IBooker& ibook,
568  edm::Run const& iRun,
569  edm::EventSetup const& iSetup) {
570  ibook.setCurrentFolder(folder_);
571  // --- histograms booking
572  meTrackEffPtTot_ = ibook.book1D("EffPtTot", "Pt of tracks associated to LV; track pt [GeV] ", 110, 0., 11.);
573  meTrackEffEtaTot_ = ibook.book1D("EffEtaTot", "Eta of tracks associated to LV; track eta ", 66, 0., 3.3);
574  meTrackMatchedTPEffPtTot_ =
575  ibook.book1D("MatchedTPEffPtTot", "Pt of tracks associated to LV matched to TP; track pt [GeV] ", 110, 0., 11.);
576  meTrackMatchedTPEffPtMtd_ = ibook.book1D(
577  "MatchedTPEffPtMtd", "Pt of tracks associated to LV matched to TP with time; track pt [GeV] ", 110, 0., 11.);
578  meTrackMatchedTPEffEtaTot_ =
579  ibook.book1D("MatchedTPEffEtaTot", "Eta of tracks associated to LV matched to TP; track eta ", 66, 0., 3.3);
580  meTrackMatchedTPEffEtaMtd_ = ibook.book1D(
581  "MatchedTPEffEtaMtd", "Eta of tracks associated to LV matched to TP with time; track eta ", 66, 0., 3.3);
582  meTrackMatchedTPResTot_ =
583  ibook.book1D("MatchedTPTrackRes",
584  "t_{rec} - t_{sim} for tracks associated to LV matched to TP; t_{rec} - t_{sim} [ns] ",
585  120,
586  -0.15,
587  0.15);
588  meTrackResTot_ = ibook.book1D("TrackRes", "t_{rec} - t_{sim} for tracks; t_{rec} - t_{sim} [ns] ", 120, -0.15, 0.15);
589  meTrackRes_[0] = ibook.book1D(
590  "TrackRes-LowMVA", "t_{rec} - t_{sim} for tracks with MVA < 0.5; t_{rec} - t_{sim} [ns] ", 100, -1., 1.);
591  meTrackRes_[1] = ibook.book1D(
592  "TrackRes-MediumMVA", "t_{rec} - t_{sim} for tracks with 0.5 < MVA < 0.8; t_{rec} - t_{sim} [ns] ", 100, -1., 1.);
593  meTrackRes_[2] = ibook.book1D(
594  "TrackRes-HighMVA", "t_{rec} - t_{sim} for tracks with MVA > 0.8; t_{rec} - t_{sim} [ns] ", 100, -1., 1.);
595  if (optionalPlots_) {
596  meTrackResMass_[0] = ibook.book1D(
597  "TrackResMass-LowMVA", "t_{rec} - t_{est} for tracks with MVA < 0.5; t_{rec} - t_{est} [ns] ", 100, -1., 1.);
598  meTrackResMass_[1] = ibook.book1D("TrackResMass-MediumMVA",
599  "t_{rec} - t_{est} for tracks with 0.5 < MVA < 0.8; t_{rec} - t_{est} [ns] ",
600  100,
601  -1.,
602  1.);
603  meTrackResMass_[2] = ibook.book1D(
604  "TrackResMass-HighMVA", "t_{rec} - t_{est} for tracks with MVA > 0.8; t_{rec} - t_{est} [ns] ", 100, -1., 1.);
605  meTrackResMassTrue_[0] = ibook.book1D(
606  "TrackResMassTrue-LowMVA", "t_{est} - t_{sim} for tracks with MVA < 0.5; t_{est} - t_{sim} [ns] ", 100, -1., 1.);
607  meTrackResMassTrue_[1] = ibook.book1D("TrackResMassTrue-MediumMVA",
608  "t_{est} - t_{sim} for tracks with 0.5 < MVA < 0.8; t_{est} - t_{sim} [ns] ",
609  100,
610  -1.,
611  1.);
612  meTrackResMassTrue_[2] = ibook.book1D("TrackResMassTrue-HighMVA",
613  "t_{est} - t_{sim} for tracks with MVA > 0.8; t_{est} - t_{sim} [ns] ",
614  100,
615  -1.,
616  1.);
617  }
618  meTrackMatchedTPPullTot_ = ibook.book1D(
619  "MatchedTPTrackPull", "Pull for tracks associated to LV matched to TP; (t_{rec}-t_{sim})/#sigma_{t}", 50, -5., 5.);
620  meTrackPullTot_ = ibook.book1D("TrackPull", "Pull for tracks; (t_{rec}-t_{sim})/#sigma_{t}", 100, -10., 10.);
621  meTrackPull_[0] =
622  ibook.book1D("TrackPull-LowMVA", "Pull for tracks with MVA < 0.5; (t_{rec}-t_{sim})/#sigma_{t}", 100, -10., 10.);
623  meTrackPull_[1] = ibook.book1D(
624  "TrackPull-MediumMVA", "Pull for tracks with 0.5 < MVA < 0.8; (t_{rec}-t_{sim})/#sigma_{t}", 100, -10., 10.);
625  meTrackPull_[2] =
626  ibook.book1D("TrackPull-HighMVA", "Pull for tracks with MVA > 0.8; (t_{rec}-t_{sim})/#sigma_{t}", 100, -10., 10.);
627  meTrackMatchedTPZposResTot_ =
628  ibook.book1D("MatchedTPTrackZposResTot",
629  "Z_{PCA} - Z_{sim} for tracks associated to LV matched to TP;Z_{PCA} - Z_{sim} [cm] ",
630  50,
631  -0.1,
632  0.1);
633  meTrackZposResTot_ =
634  ibook.book1D("TrackZposResTot", "Z_{PCA} - Z_{sim} for tracks;Z_{PCA} - Z_{sim} [cm] ", 50, -0.5, 0.5);
635  meTrackZposRes_[0] = ibook.book1D(
636  "TrackZposRes-LowMVA", "Z_{PCA} - Z_{sim} for tracks with MVA < 0.5;Z_{PCA} - Z_{sim} [cm] ", 50, -0.5, 0.5);
637  meTrackZposRes_[1] = ibook.book1D("TrackZposRes-MediumMVA",
638  "Z_{PCA} - Z_{sim} for tracks with 0.5 < MVA < 0.8 ;Z_{PCA} - Z_{sim} [cm] ",
639  50,
640  -0.5,
641  0.5);
642  meTrackZposRes_[2] = ibook.book1D(
643  "TrackZposRes-HighMVA", "Z_{PCA} - Z_{sim} for tracks with MVA > 0.8 ;Z_{PCA} - Z_{sim} [cm] ", 50, -0.5, 0.5);
644  meTrack3DposRes_[0] =
645  ibook.book1D("Track3DposRes-LowMVA",
646  "3dPos_{PCA} - 3dPos_{sim} for tracks with MVA < 0.5 ;3dPos_{PCA} - 3dPos_{sim} [cm] ",
647  50,
648  -0.5,
649  0.5);
650  meTrack3DposRes_[1] =
651  ibook.book1D("Track3DposRes-MediumMVA",
652  "3dPos_{PCA} - 3dPos_{sim} for tracks with 0.5 < MVA < 0.8 ;3dPos_{PCA} - 3dPos_{sim} [cm] ",
653  50,
654  -0.5,
655  0.5);
656  meTrack3DposRes_[2] =
657  ibook.book1D("Track3DposRes-HighMVA",
658  "3dPos_{PCA} - 3dPos_{sim} for tracks with MVA > 0.8;3dPos_{PCA} - 3dPos_{sim} [cm] ",
659  50,
660  -0.5,
661  0.5);
662  meTimeRes_ = ibook.book1D("TimeRes", "t_{rec} - t_{sim} ;t_{rec} - t_{sim} [ns] ", 40, -0.2, 0.2);
663  meTimePull_ = ibook.book1D("TimePull", "Pull; t_{rec} - t_{sim}/#sigma_{t rec}", 100, -10., 10.);
664  meTimeSignalRes_ =
665  ibook.book1D("TimeSignalRes", "t_{rec} - t_{sim} for signal ;t_{rec} - t_{sim} [ns] ", 40, -0.2, 0.2);
666  meTimeSignalPull_ =
667  ibook.book1D("TimeSignalPull", "Pull for signal; t_{rec} - t_{sim}/#sigma_{t rec}", 100, -10., 10.);
668  mePUvsRealV_ =
669  ibook.bookProfile("PUvsReal", "#PU vertices vs #real matched vertices;#PU;#real ", 100, 0, 300, 100, 0, 200);
670  mePUvsFakeV_ =
671  ibook.bookProfile("PUvsFake", "#PU vertices vs #fake matched vertices;#PU;#fake ", 100, 0, 300, 100, 0, 20);
672  mePUvsOtherFakeV_ = ibook.bookProfile(
673  "PUvsOtherFake", "#PU vertices vs #other fake matched vertices;#PU;#other fake ", 100, 0, 300, 100, 0, 20);
674  mePUvsSplitV_ =
675  ibook.bookProfile("PUvsSplit", "#PU vertices vs #split matched vertices;#PU;#split ", 100, 0, 300, 100, 0, 20);
676  meMatchQual_ = ibook.book1D("MatchQuality", "RECO-SIM vertex match quality; ", 8, 0, 8.);
677  meDeltaZrealreal_ = ibook.book1D("DeltaZrealreal", "#Delta Z real-real; |#Delta Z (r-r)| [cm]", 100, 0, 0.5);
678  meDeltaZfakefake_ = ibook.book1D("DeltaZfakefake", "#Delta Z fake-fake; |#Delta Z (f-f)| [cm]", 100, 0, 0.5);
679  meDeltaZfakereal_ = ibook.book1D("DeltaZfakereal", "#Delta Z fake-real; |#Delta Z (f-r)| [cm]", 100, 0, 0.5);
680  meDeltaTrealreal_ = ibook.book1D("DeltaTrealreal", "#Delta T real-real; |#Delta T (r-r)| [sigma]", 60, 0., 30.);
681  meDeltaTfakefake_ = ibook.book1D("DeltaTfakefake", "#Delta T fake-fake; |#Delta T (f-f)| [sigma]", 60, 0., 30.);
682  meDeltaTfakereal_ = ibook.book1D("DeltaTfakereal", "#Delta T fake-real; |#Delta T (f-r)| [sigma]", 60, 0., 30.);
683  if (optionalPlots_) {
684  meRecoPosInSimCollection_ = ibook.book1D(
685  "RecoPosInSimCollection", "Sim signal vertex index associated to Reco signal vertex; Sim PV index", 200, 0, 200);
686  meRecoPosInRecoOrigCollection_ =
687  ibook.book1D("RecoPosInRecoOrigCollection", "Reco signal index in OrigCollection; Reco index", 200, 0, 200);
688  meSimPosInSimOrigCollection_ =
689  ibook.book1D("SimPosInSimOrigCollection", "Sim signal index in OrigCollection; Sim index", 200, 0, 200);
690  }
691  meRecoPVPosSignal_ =
692  ibook.book1D("RecoPVPosSignal", "Position in reco collection of PV associated to sim signal", 20, 0, 20);
693  meRecoPVPosSignalNotHighestPt_ =
694  ibook.book1D("RecoPVPosSignalNotHighestPt",
695  "Position in reco collection of PV associated to sim signal not highest Pt",
696  20,
697  0,
698  20);
699  meRecoVtxVsLineDensity_ =
700  ibook.book1D("RecoVtxVsLineDensity", "#Reco vertices/mm/event; line density [#vtx/mm/event]", 160, 0., 4.);
701  meRecVerNumber_ = ibook.book1D("RecVerNumber", "RECO Vertex Number: Number of vertices", 50, 0, 250);
702  meRecPVZ_ = ibook.book1D("recPVZ", "Weighted #Rec vertices/mm", 400, -20., 20.);
703  meRecPVT_ = ibook.book1D("recPVT", "#Rec vertices/10 ps", 200, -1., 1.);
704  meSimPVZ_ = ibook.book1D("simPVZ", "Weighted #Sim vertices/mm", 400, -20., 20.);
705 
706  meVtxTrackMult_ = ibook.book1D("VtxTrackMult", "Log10(Vertex track multiplicity)", 80, 0.5, 2.5);
707  meVtxTrackW_ = ibook.book1D("VtxTrackW", "Vertex track weight (all)", 50, 0., 1.);
708  meVtxTrackWnt_ = ibook.book1D("VtxTrackWnt", "Vertex track Wnt", 50, 0., 1.);
709  meVtxTrackRecLVMult_ =
710  ibook.book1D("VtxTrackRecLVMult", "Log10(Vertex track multiplicity) for matched LV", 80, 0.5, 2.5);
711  meVtxTrackRecLVW_ = ibook.book1D("VtxTrackRecLVW", "Vertex track weight for matched LV (all)", 50, 0., 1.);
712  meVtxTrackRecLVWnt_ = ibook.book1D("VtxTrackRecLVWnt", "Vertex track Wnt for matched LV", 50, 0., 1.);
713 
714  mePUTrackRelMult_ = ibook.book1D(
715  "PUTrackRelMult", "Relative multiplicity of PU tracks for matched vertices; #PUTrks/#Trks", 50, 0., 1.);
716  meFakeTrackRelMult_ = ibook.book1D(
717  "FakeTrackRelMult", "Relative multiplicity of fake tracks for matched vertices; #fakeTrks/#Trks", 50, 0., 1.);
718  mePUTrackRelSumWnt_ =
719  ibook.book1D("PUTrackRelSumWnt",
720  "Relative Sum of wnt of PU tracks for matched vertices; PUSumW*min(Pt, 1.)/SumW*min(Pt, 1.)",
721  50,
722  0.,
723  1.);
724  mePUTrackRelSumWos_ = ibook.book1D(
725  "PUTrackRelSumWos", "Relative Sum of wos of PU tracks for matched vertices; PUSumWos/SumWos", 50, 0., 1.);
726  meSecTrackRelSumWos_ =
727  ibook.book1D("SecTrackRelSumWos",
728  "Relative Sum of wos of tracks from secondary vtx for matched vertices; SecSumWos/SumWos",
729  50,
730  0.,
731  1.);
732  meFakeTrackRelSumWos_ = ibook.book1D(
733  "FakeTrackRelSumWos", "Relative Sum of wos of fake tracks for matched vertices; FakeSumWos/SumWos", 50, 0., 1.);
734  mePUTrackRelSumPt_ = ibook.book1D(
735  "PUTrackRelSumPt", "Relative Sum of Pt of PU tracks for matched vertices; PUSumPt/SumPt", 50, 0., 1.);
736  mePUTrackRelSumPt2_ = ibook.book1D(
737  "PUTrackRelSumPt2", "Relative Sum of Pt2 for PU tracks for matched vertices; PUSumPt2/SumPt2", 50, 0., 1.);
738  mePUTrackRecLVRelMult_ = ibook.book1D(
739  "PUTrackRecLVRelMult", "Relative multiplicity of PU tracks for matched LV; #PUTrks/#Trks", 50, 0., 1.);
740  meFakeTrackRecLVRelMult_ = ibook.book1D(
741  "FakeTrackRecLVRelMult", "Relative multiplicity of fake tracks for matched LV; #FakeTrks/#Trks", 50, 0., 1.);
742  mePUTrackRecLVRelSumWnt_ =
743  ibook.book1D("PUTrackRecLVRelSumWnt",
744  "Relative Sum of Wnt of PU tracks for matched LV; PUSumW*min(Pt, 1.)/SumW*min(Pt, 1.)",
745  50,
746  0.,
747  1.);
748  mePUTrackRecLVRelSumWos_ = ibook.book1D(
749  "PUTrackRecLVRelSumWos", "Relative Sum of Wos of PU tracks for matched LV; PUSumWos/SumWos", 50, 0., 1.);
750  meSecTrackRecLVRelSumWos_ =
751  ibook.book1D("SecTrackRecLVRelSumWos",
752  "Relative Sum of wos of tracks from secondary vtx for matched LV; SecSumWos/SumWos",
753  50,
754  0.,
755  1.);
756  meFakeTrackRecLVRelSumWos_ = ibook.book1D(
757  "FakeTrackRecLVRelSumWos", "Relative Sum of wos of fake tracks for matched LV; FakeSumWos/SumWos", 50, 0., 1.);
758  mePUTrackRecLVRelSumPt_ =
759  ibook.book1D("PUTrackRecLVRelSumPt", "Relative Sum of Pt of PU tracks for matched LV; PUSumPt/SumPt", 50, 0., 1.);
760  mePUTrackRecLVRelSumPt2_ = ibook.book1D(
761  "PUTrackRecLVRelSumPt2", "Relative Sum of Pt2 of PU tracks for matched LV; PUSumPt2/SumPt2", 50, 0., 1.);
762 
763  if (optionalPlots_) {
764  mePUTrackMult_ = ibook.book1D("PUTrackMult", "Number of PU tracks for matched vertices; #PUTrks", 50, 0., 100.);
765  mePUTrackWnt_ = ibook.book1D("PUTrackWnt", "Wnt of PU tracks for matched vertices; PUTrkW*min(Pt, 1.)", 50, 0., 1.);
766  mePUTrackSumWnt_ = ibook.book1D(
767  "PUTrackSumWnt", "Sum of wnt of PU tracks for matched vertices; log10(PUSumW*min(Pt, 1.))", 50, -2., 3.);
768  mePUTrackSumWos_ =
769  ibook.book1D("PUTrackSumWos", "Sum of wos of PU tracks for matched vertices; log10(PUSumWos)", 50, -1., 7.);
770  meSecTrackSumWos_ = ibook.book1D(
771  "SecTrackSumWos", "Sum of wos of tracks from secondary vtx for matched vertices; log10(SecSumWos)", 50, -1., 7.);
772  mePUTrackSumPt_ =
773  ibook.book1D("PUTrackSumPt", "Sum of Pt of PU tracks for matched vertices; log10(PUSumPt)", 50, -2., 3.);
774  mePUTrackSumPt2_ =
775  ibook.book1D("PUTrackSumPt2", "Sum of Pt2 of PU tracks for matched vertices; log10(PUSumPt2)", 50, -3., 3.);
776 
777  mePUTrackRelMultvsMult_ =
778  ibook.bookProfile("PUTrackRelMultvsMult",
779  "Relative PU multiplicity vs Number of tracks for matched vertices; #Trks; #PUTrks/#Trks",
780  50,
781  0.,
782  200.,
783  0.,
784  1.,
785  "s");
786  meFakeTrackRelMultvsMult_ = ibook.bookProfile(
787  "FakeTrackRelMultvsMult",
788  "Relative multiplicity of fake tracks vs Number of tracks for matched vertices; #Trks; #FakeTrks/#Trks",
789  50,
790  0.,
791  200.,
792  0.,
793  1.,
794  "s");
795  mePUTrackRelSumWntvsSumWnt_ =
796  ibook.bookProfile("PUTrackRelSumWntvsSumWnt",
797  "Relative PU Sum of Wnt vs Sum of Wnt of tracks for matched vertices; log10(SumW*min(Pt, "
798  "1.)); PUSumW*min(Pt, 1.)/SumW*min(Pt, 1.)",
799  50,
800  0.,
801  2.5,
802  0.,
803  1.,
804  "s");
805  mePUTrackRelSumWosvsSumWos_ = ibook.bookProfile(
806  "PUTrackRelSumWosvsSumWos",
807  "Relative PU Sum of Wos vs Sum of Wos of tracks for matched vertices; log10(SumWos); PUSumWos/SumWos",
808  50,
809  2.5,
810  7.,
811  0.,
812  1.,
813  "s");
814  meSecTrackRelSumWosvsSumWos_ = ibook.bookProfile("SecTrackRelSumWosvsSumWos",
815  "Relative Sum of Wos of tracks from secondary vtx vs Sum of Wos "
816  "of tracks for matched vertices; log10(SumWos); SecSumWos/SumWos",
817  50,
818  2.,
819  7.,
820  0.,
821  1.,
822  "s");
823  meFakeTrackRelSumWosvsSumWos_ = ibook.bookProfile("FakeTrackRelSumWosvsSumWos",
824  "Relative Sum of Wos of fake tracks vs Sum of Wos of tracks for "
825  "matched vertices; log10(SumWos); FakeSumWos/SumWos",
826  50,
827  2.5,
828  7.5,
829  0.,
830  1.,
831  "s");
832  mePUTrackRelSumPtvsSumPt_ = ibook.bookProfile(
833  "PUTrackRelSumPtvsSumPt",
834  "Relative PU Sum of Pt vs Sum of Pt of tracks for matched vertices; log10(SumPt); PUSumPt/SumPt",
835  50,
836  0.,
837  3.,
838  0.,
839  1.,
840  "s");
841  mePUTrackRelSumPt2vsSumPt2_ = ibook.bookProfile(
842  "PUTrackRelSumPt2vsSumPt2",
843  "Relative PU Sum of Pt2 vs Sum of Pt2 of tracks for matched vertices; log10(SumPt2); PUSumPt2/SumPt2",
844  50,
845  0.,
846  4.,
847  0.,
848  1.,
849  "s");
850 
851  mePUTrackRecLVMult_ = ibook.book1D("PUTrackRecLVMult", "Number of PU tracks for matched LV; #PUTrks", 50, 0., 100.);
852  mePUTrackRecLVWnt_ =
853  ibook.book1D("PUTrackRecLVWnt", "Wnt of PU tracks for matched LV; PUTrkW*min(Pt, 1.)", 50, 0., 1.);
854  mePUTrackRecLVSumWnt_ = ibook.book1D(
855  "PUTrackRecLVSumWnt", "Sum of wnt of PU tracks for matched LV; log10(PUSumW*min(Pt, 1.))", 50, -2., 3.);
856  mePUTrackRecLVSumWos_ =
857  ibook.book1D("PUTrackRecLVSumWos", "Sum of wos of PU tracks for matched LV; log10(PUSumWos)", 50, -1., 7.);
858  meSecTrackRecLVSumWos_ = ibook.book1D(
859  "SecTrackRecLVSumWos", "Sum of wos of tracks from secondary vtx for matched LV; log10(SecSumWos)", 50, -1., 7.);
860  mePUTrackRecLVSumPt_ =
861  ibook.book1D("PUTrackRecLVSumPt", "Sum of Pt of PU tracks for matched LV; log10(PUSumPt)", 50, -2., 3.);
862  mePUTrackRecLVSumPt2_ =
863  ibook.book1D("PUTrackRecLVSumPt2", "Sum of Pt2 of PU tracks for matched LV; log10(PUSumPt2)", 50, -3., 3.);
864 
865  mePUTrackRecLVRelMultvsMult_ =
866  ibook.bookProfile("PUTrackRecLVRelMultvsMult",
867  "Relative PU multiplicity vs Number of tracks for matched LV; #Trks; #PUTrks/#Trks",
868  50,
869  0.,
870  200.,
871  0.,
872  1.,
873  "s");
874  meFakeTrackRecLVRelMultvsMult_ = ibook.bookProfile(
875  "FakeTrackRecLVRelMultvsMult",
876  "Relative multiplicity of fake tracks vs Number of tracks for matched LV; #Trks; #FakeTrks/#Trks",
877  50,
878  0.,
879  200.,
880  0.,
881  1.,
882  "s");
883  mePUTrackRecLVRelSumWntvsSumWnt_ =
884  ibook.bookProfile("PUTrackRecLVRelSumWntvsSumWnt",
885  "Relative PU Sum of Wnt vs Sum of Wnt of tracks for matched LV; log10(SumW*min(Pt, 1.)); "
886  "PUSumW*min(Pt, 1.)/SumW*min(Pt, 1.)",
887  50,
888  1.,
889  2.3,
890  0.,
891  1.,
892  "s");
893  mePUTrackRecLVRelSumWosvsSumWos_ = ibook.bookProfile(
894  "PUTrackRecLVRelSumWosvsSumWos",
895  "Relative PU Sum of Wos vs Sum of Wos of tracks for matched vertices; log10(SumWos); PUSumWos/SumWos",
896  50,
897  5.5,
898  7.,
899  0.,
900  1.,
901  "s");
902  meSecTrackRecLVRelSumWosvsSumWos_ =
903  ibook.bookProfile("SecTrackRecLVRelSumWosvsSumWos",
904  "Relative Sum of Wos of tracks from secondary vtx vs Sum of Wos of tracks for matched "
905  "vertices; log10(SumWos); SecSumWos/SumWos",
906  50,
907  5.5,
908  7.,
909  0.,
910  1.,
911  "s");
912  meFakeTrackRecLVRelSumWosvsSumWos_ = ibook.bookProfile("FakeTrackRecLVRelSumWosvsSumWos",
913  "Relative Sum of Wos of fake tracks vs Sum of Wos of tracks "
914  "for matched vertices; log10(SumWos); FakeSumWos/SumWos",
915  50,
916  5.5,
917  7.,
918  0.,
919  1.,
920  "s");
921  mePUTrackRecLVRelSumPtvsSumPt_ =
922  ibook.bookProfile("PUTrackRecLVRelSumPtvsSumPt",
923  "Relative PU Sum of Pt vs Sum of Pt of tracks for matched LV; log10(SumPt); PUSumPt/SumPt",
924  50,
925  1.4,
926  3.,
927  0.,
928  1.,
929  "s");
930  mePUTrackRecLVRelSumPt2vsSumPt2_ =
931  ibook.bookProfile("PUTrackRecLVRelSumPt2vsSumPt2",
932  "Relative PU Sum of Pt2 vs Sum of tracks for matched LV; log10(SumPt2); PUSumPt2/SumPt2",
933  50,
934  2.,
935  4.,
936  0.,
937  1.,
938  "s");
939  }
940 
941  meJetsPURelMult_ =
942  ibook.book1D("JetsPURelMult",
943  "Relative contribution of PU to jet multiplicity for matched vertices; (#Jets-#JetsNoPU)/#Jets",
944  50,
945  0.,
946  1.);
947  meJetsPURelHt_ =
948  ibook.book1D("JetsPURelHt",
949  "Relative contribution of PU to scalar sum of Et of jets for matched vertices; (Ht-HtNoPU)/HT",
950  50,
951  0.,
952  1.);
953  meJetsPURelSumPt2_ =
954  ibook.book1D("JetsPURelSumPt2",
955  "Relative contribution of PU to sum of Pt2 of jets for matched vertices; (SumPt2-SumPt2NoPU)/SumPt2",
956  50,
957  0.,
958  1.);
959  meJetsFakeRelSumPt2_ = ibook.book1D(
960  "JetsFakeRelSumPt2",
961  "Relative contribution of fake tracks to sum of Pt2 of jets for matched vertices; (SumPt2-SumPt2NoFake)/SumPt2",
962  50,
963  0.,
964  1.);
965  meJetsPURelMetPt_ =
966  ibook.book1D("JetsPURelMetPt",
967  "Relative contribution of PU to Missing Transverse Energy for matched vertices; (Met-MetNoPU)/Met",
968  50,
969  -1.,
970  1.);
971  meJetsPURelSumPz_ =
972  ibook.book1D("JetsPURelSumPz",
973  "Relative contribution of PU to sum of Pz of jets for matched vertices; (SumPz-SumPzNoPU)/SumPz",
974  50,
975  -1.,
976  1.);
977 
978  meJetsRecLVPURelMult_ =
979  ibook.book1D("JetsRecLVPURelMult",
980  "Relative contribution of PU to jet multiplicity for matched LV; (#Jets-#JetsNoPU)/#Jets",
981  50,
982  0.,
983  1.);
984  meJetsRecLVPURelHt_ =
985  ibook.book1D("JetsRecLVPURelHt",
986  "Relative contribution of PU to scalar sum of Et of jets for matched LV; (Ht-HtNoPU)/HT",
987  50,
988  0.,
989  1.);
990  meJetsRecLVPURelSumPt2_ =
991  ibook.book1D("JetsRecLVPURelSumPt2",
992  "Relative contribution of PU to sum of Pt2 of jets for matched LV; (SumPt2-SumPt2NoPU)/SumPt2",
993  50,
994  0.,
995  1.);
996  meJetsRecLVFakeRelSumPt2_ = ibook.book1D(
997  "JetsRecLVFakeRelSumPt2",
998  "Relative contribution of fake tracks to sum of Pt2 of jets for matched LV; (SumPt2-SumPt2NoFake)/SumPt2",
999  50,
1000  0.,
1001  1.);
1002  meJetsRecLVPURelMetPt_ =
1003  ibook.book1D("JetsRecLVPURelMetPt",
1004  "Relative contribution of PU to Missing Transverse Energy for matched LV; (Met-MetNoPU)/Met",
1005  50,
1006  -1.,
1007  1.);
1008  meJetsRecLVPURelSumPz_ =
1009  ibook.book1D("JetsRecLVPURelSumPz",
1010  "Relative contribution of PU to sum of Pz of jets for matched LV; (SumPz-SumPzNoPU)/SumPz",
1011  50,
1012  -1.,
1013  1.);
1014 
1015  if (optionalPlots_) {
1016  meJetsPUMult_ = ibook.book1D(
1017  "JetsPUMult", "Contribution of PU to jet multiplicity for matched vertices; #Jets-#JetsNoPU", 50, 0., 100.);
1018  meJetsPUHt_ = ibook.book1D("JetsPUHt",
1019  "Contribution of PU to scalar sum of Et of jets for matched vertices; log10(Ht-HtNoPU)",
1020  50,
1021  -2.,
1022  3.);
1023  meJetsPUSumPt2_ =
1024  ibook.book1D("JetsPUSumPt2",
1025  "Contribution of PU to sum of Pt2 of jets for matched vertices; log10(sumPt2-SumPt2NoPU)",
1026  50,
1027  -3.,
1028  3.);
1029  meJetsPUMetPt_ =
1030  ibook.book1D("JetsPUMetPt",
1031  "Contribution of PU to Missing Transverse Energy for matched vertices; log10(Met-MetNoPU)",
1032  50,
1033  -3.,
1034  2.);
1035  meJetsPUSumPz_ =
1036  ibook.book1D("JetsPUSumPz",
1037  "Contribution of PU to sum of Pz of jets for matched vertices; log10(abs(SumPz-SumPzNoPU))",
1038  50,
1039  -4.,
1040  3.);
1041 
1042  meJetsPURelMultvsMult_ = ibook.bookProfile("JetsPURelMultvsMult",
1043  "Relative contribution of PU to jet multiplicity vs number of jets for "
1044  "matched vertices; #Jets; (#Jets-#JetsNoPU)/#Jets",
1045  50,
1046  0.,
1047  120.,
1048  0.,
1049  1.,
1050  "s");
1051  meJetsPURelHtvsHt_ = ibook.bookProfile("JetsPURelHtvsHt",
1052  "Relative contribution of PU to scalar sum of Et of jets vs scalar sum of "
1053  "Et for matched vertices; log10(Ht); (Ht-HtNoPU)/HT",
1054  50,
1055  0.,
1056  3.,
1057  0.,
1058  1.,
1059  "s");
1060  meJetsPURelSumPt2vsSumPt2_ = ibook.bookProfile("JetsPURelSumPt2vsSumPt2",
1061  "Relative contribution of PU to sum of Pt2 of jets vs sum of Pt2 "
1062  "for matched vertices; log10(SumPt2); (SumPt2-SumPt2NoPU)/SumPt2",
1063  50,
1064  -1.,
1065  4.,
1066  0.,
1067  1.,
1068  "s");
1069  meJetsFakeRelSumPt2vsSumPt2_ =
1070  ibook.bookProfile("JetsFakeRelSumPt2vsSumPt2",
1071  "Relative contribution of fake tracks to sum of Pt2 of jets vs sum of Pt2 for matched "
1072  "vertices; log10(SumPt2); (SumPt2-SumPt2NoFake)/SumPt2",
1073  50,
1074  -1.,
1075  4.,
1076  0.,
1077  1.,
1078  "s");
1079  meJetsPURelMetPtvsMetPt_ = ibook.bookProfile("JetsPURelMetPtvsMetPt",
1080  "Relative contribution of PU to Missing Transverse Energy vs MET for "
1081  "matched vertices; log10(Met); (Met-MetNoPU)/Met",
1082  50,
1083  -1.,
1084  2.,
1085  -1.,
1086  1.,
1087  "s");
1088  meJetsPURelSumPzvsSumPz_ = ibook.bookProfile("JetsPURelSumPzvsSumPz",
1089  "Relative contribution of PU to sum of Pz of jets vs Sum of Pz for "
1090  "matched vertices; log10(abs SumPz); (SumPz-SumPzNoPU)/SumPz",
1091  50,
1092  -2.,
1093  3.,
1094  -1.,
1095  1.,
1096  "s");
1097 
1098  meJetsRecLVPUMult_ = ibook.book1D(
1099  "JetsRecLVPUMult", "Contribution of PU to jet multiplicity for matched LV; #Jets-#JetsNoPU", 50, 0., 100.);
1100  meJetsRecLVPUHt_ = ibook.book1D(
1101  "JetsRecLVPUHt", "Contribution of PU to scalar sum of Et of jets for matched LV; log10(Ht-HtNoPU)", 50, -2., 3.);
1102  meJetsRecLVPUSumPt2_ =
1103  ibook.book1D("JetsRecLVPUSumPt2",
1104  "Contribution of PU to sum of Pt2 of jets for matched LV; log10(sumPt2-SumPt2NoPU)",
1105  50,
1106  -3.,
1107  3.);
1108  meJetsRecLVPUMetPt_ =
1109  ibook.book1D("JetsRecLVPUMetPt",
1110  "Contribution of PU to Missing Transverse Energy for matched LV; log10(Met-MetNoPU)",
1111  50,
1112  -3.,
1113  2.);
1114  meJetsRecLVPUSumPz_ =
1115  ibook.book1D("JetsRecLVPUSumPz",
1116  "Contribution of PU to sum of Pz of jets for matched LV; log10(abs(SumPz-SumPzNoPU))",
1117  50,
1118  -4.,
1119  3.);
1120 
1121  meJetsRecLVPURelMultvsMult_ = ibook.bookProfile("JetsRecLVPURelMultvsMult",
1122  "Relative contribution of PU to jet multiplicity vs number of jets "
1123  "for matched vertices; #Jets; (#Jets-#JetsNoPU)/#Jets",
1124  50,
1125  0.,
1126  120.,
1127  0.,
1128  1.,
1129  "s");
1130  meJetsRecLVPURelHtvsHt_ = ibook.bookProfile("JetsRecLVPURelHtvsHt",
1131  "Relative contribution of PU to scalar sum of Et of jets vs scalar sum "
1132  "of Et for matched vertices; log10(Ht); (Ht-HtNoPU)/HT",
1133  50,
1134  1.5,
1135  3.,
1136  0.,
1137  1.,
1138  "s");
1139  meJetsRecLVPURelSumPt2vsSumPt2_ =
1140  ibook.bookProfile("JetsRecLVPURelSumPt2vsSumPt2",
1141  "Relative contribution of PU to sum of Pt2 of jets vs sum of Pt2 for matched vertices; "
1142  "log10(SumPt2); (SumPt2-SumPt2NoPU)/SumPt2",
1143  50,
1144  2.,
1145  5.,
1146  0.,
1147  1.,
1148  "s");
1149  meJetsRecLVFakeRelSumPt2vsSumPt2_ =
1150  ibook.bookProfile("JetsRecLVFakeRelSumPt2vsSumPt2",
1151  "Relative contribution of fake tracks to sum of Pt2 of jets vs sum of Pt2 for matched "
1152  "vertices; log10(SumPt2); (SumPt2-SumPt2NoFake)/SumPt2",
1153  50,
1154  2.,
1155  5.,
1156  0.,
1157  1.,
1158  "s");
1159  meJetsRecLVPURelMetPtvsMetPt_ = ibook.bookProfile("JetsRecLVPURelMetPtvsMetPt",
1160  "Relative contribution of PU to Missing Transverse Energy vs MET "
1161  "for matched vertices; log10(Met); (Met-MetNoPU)/Met",
1162  50,
1163  0.,
1164  2.5,
1165  -1.,
1166  1.,
1167  "s");
1168  meJetsRecLVPURelSumPzvsSumPz_ = ibook.bookProfile("JetsRecLVPURelSumPzvsSumPz",
1169  "Relative contribution of PU to sum of Pz of jets vs Sum of Pz "
1170  "for matched vertices; log10(abs SumPz); (SumPz-SumPzNoPU)/SumPz",
1171  50,
1172  0.5,
1173  3.5,
1174  -1.,
1175  1.,
1176  "s");
1177  }
1178 
1179  // some tests
1180  meTrackResLowPTot_ = ibook.book1D(
1181  "TrackResLowP", "t_{rec} - t_{sim} for tracks with p < 2 GeV; t_{rec} - t_{sim} [ns] ", 70, -0.15, 0.15);
1182  meTrackResLowP_[0] =
1183  ibook.book1D("TrackResLowP-LowMVA",
1184  "t_{rec} - t_{sim} for tracks with MVA < 0.5 and p < 2 GeV; t_{rec} - t_{sim} [ns] ",
1185  100,
1186  -1.,
1187  1.);
1188  meTrackResLowP_[1] =
1189  ibook.book1D("TrackResLowP-MediumMVA",
1190  "t_{rec} - t_{sim} for tracks with 0.5 < MVA < 0.8 and p < 2 GeV; t_{rec} - t_{sim} [ns] ",
1191  100,
1192  -1.,
1193  1.);
1194  meTrackResLowP_[2] =
1195  ibook.book1D("TrackResLowP-HighMVA",
1196  "t_{rec} - t_{sim} for tracks with MVA > 0.8 and p < 2 GeV; t_{rec} - t_{sim} [ns] ",
1197  100,
1198  -1.,
1199  1.);
1200  meTrackResHighPTot_ = ibook.book1D(
1201  "TrackResHighP", "t_{rec} - t_{sim} for tracks with p > 2 GeV; t_{rec} - t_{sim} [ns] ", 70, -0.15, 0.15);
1202  meTrackResHighP_[0] =
1203  ibook.book1D("TrackResHighP-LowMVA",
1204  "t_{rec} - t_{sim} for tracks with MVA < 0.5 and p > 2 GeV; t_{rec} - t_{sim} [ns] ",
1205  100,
1206  -1.,
1207  1.);
1208  meTrackResHighP_[1] =
1209  ibook.book1D("TrackResHighP-MediumMVA",
1210  "t_{rec} - t_{sim} for tracks with 0.5 < MVA < 0.8 and p > 2 GeV; t_{rec} - t_{sim} [ns] ",
1211  100,
1212  -1.,
1213  1.);
1214  meTrackResHighP_[2] =
1215  ibook.book1D("TrackResHighP-HighMVA",
1216  "t_{rec} - t_{sim} for tracks with MVA > 0.8 and p > 2 GeV; t_{rec} - t_{sim} [ns] ",
1217  100,
1218  -1.,
1219  1.);
1220  meTrackPullLowPTot_ =
1221  ibook.book1D("TrackPullLowP", "Pull for tracks with p < 2 GeV; (t_{rec}-t_{sim})/#sigma_{t}", 100, -10., 10.);
1222  meTrackPullLowP_[0] = ibook.book1D("TrackPullLowP-LowMVA",
1223  "Pull for tracks with MVA < 0.5 and p < 2 GeV; (t_{rec}-t_{sim})/#sigma_{t}",
1224  100,
1225  -10.,
1226  10.);
1227  meTrackPullLowP_[1] = ibook.book1D("TrackPullLowP-MediumMVA",
1228  "Pull for tracks with 0.5 < MVA < 0.8 and p < 2 GeV; (t_{rec}-t_{sim})/#sigma_{t}",
1229  100,
1230  -10.,
1231  10.);
1232  meTrackPullLowP_[2] = ibook.book1D("TrackPullLowP-HighMVA",
1233  "Pull for tracks with MVA > 0.8 and p < 2 GeV; (t_{rec}-t_{sim})/#sigma_{t}",
1234  100,
1235  -10.,
1236  10.);
1237  meTrackPullHighPTot_ =
1238  ibook.book1D("TrackPullHighP", "Pull for tracks with p > 2 GeV; (t_{rec}-t_{sim})/#sigma_{t}", 100, -10., 10.);
1239  meTrackPullHighP_[0] = ibook.book1D("TrackPullHighP-LowMVA",
1240  "Pull for tracks with MVA < 0.5 and p > 2 GeV; (t_{rec}-t_{sim})/#sigma_{t}",
1241  100,
1242  -10.,
1243  10.);
1244  meTrackPullHighP_[1] =
1245  ibook.book1D("TrackPullHighP-MediumMVA",
1246  "Pull for tracks with 0.5 < MVA < 0.8 and p > 2 GeV; (t_{rec}-t_{sim})/#sigma_{t}",
1247  100,
1248  -10.,
1249  10.);
1250  meTrackPullHighP_[2] = ibook.book1D("TrackPullHighP-HighMVA",
1251  "Pull for tracks with MVA > 0.8 and p > 2 GeV; (t_{rec}-t_{sim})/#sigma_{t}",
1252  100,
1253  -10.,
1254  10.);
1255  if (optionalPlots_) {
1256  meTrackResMassProtons_[0] =
1257  ibook.book1D("TrackResMass-Protons-LowMVA",
1258  "t_{rec} - t_{est} for proton tracks with MVA < 0.5; t_{rec} - t_{est} [ns] ",
1259  100,
1260  -1.,
1261  1.);
1262  meTrackResMassProtons_[1] =
1263  ibook.book1D("TrackResMass-Protons-MediumMVA",
1264  "t_{rec} - t_{est} for proton tracks with 0.5 < MVA < 0.8; t_{rec} - t_{est} [ns] ",
1265  100,
1266  -1.,
1267  1.);
1268  meTrackResMassProtons_[2] =
1269  ibook.book1D("TrackResMass-Protons-HighMVA",
1270  "t_{rec} - t_{est} for proton tracks with MVA > 0.8; t_{rec} - t_{est} [ns] ",
1271  100,
1272  -1.,
1273  1.);
1274  meTrackResMassTrueProtons_[0] =
1275  ibook.book1D("TrackResMassTrue-Protons-LowMVA",
1276  "t_{est} - t_{sim} for proton tracks with MVA < 0.5; t_{est} - t_{sim} [ns] ",
1277  100,
1278  -1.,
1279  1.);
1280  meTrackResMassTrueProtons_[1] =
1281  ibook.book1D("TrackResMassTrue-Protons-MediumMVA",
1282  "t_{est} - t_{sim} for proton tracks with 0.5 < MVA < 0.8; t_{est} - t_{sim} [ns] ",
1283  100,
1284  -1.,
1285  1.);
1286  meTrackResMassTrueProtons_[2] =
1287  ibook.book1D("TrackResMassTrue-Protons-HighMVA",
1288  "t_{est} - t_{sim} for proton tracks with MVA > 0.8; t_{est} - t_{sim} [ns] ",
1289  100,
1290  -1.,
1291  1.);
1292 
1293  meTrackResMassPions_[0] = ibook.book1D("TrackResMass-Pions-LowMVA",
1294  "t_{rec} - t_{est} for pion tracks with MVA < 0.5; t_{rec} - t_{est} [ns] ",
1295  100,
1296  -1.,
1297  1.);
1298  meTrackResMassPions_[1] =
1299  ibook.book1D("TrackResMass-Pions-MediumMVA",
1300  "t_{rec} - t_{est} for pion tracks with 0.5 < MVA < 0.8; t_{rec} - t_{est} [ns] ",
1301  100,
1302  -1.,
1303  1.);
1304  meTrackResMassPions_[2] = ibook.book1D("TrackResMass-Pions-HighMVA",
1305  "t_{rec} - t_{est} for pion tracks with MVA > 0.8; t_{rec} - t_{est} [ns] ",
1306  100,
1307  -1.,
1308  1.);
1309  meTrackResMassTruePions_[0] =
1310  ibook.book1D("TrackResMassTrue-Pions-LowMVA",
1311  "t_{est} - t_{sim} for pion tracks with MVA < 0.5; t_{est} - t_{sim} [ns] ",
1312  100,
1313  -1.,
1314  1.);
1315  meTrackResMassTruePions_[1] =
1316  ibook.book1D("TrackResMassTrue-Pions-MediumMVA",
1317  "t_{est} - t_{sim} for pion tracks with 0.5 < MVA < 0.8; t_{est} - t_{sim} [ns] ",
1318  100,
1319  -1.,
1320  1.);
1321  meTrackResMassTruePions_[2] =
1322  ibook.book1D("TrackResMassTrue-Pions-HighMVA",
1323  "t_{est} - t_{sim} for pion tracks with MVA > 0.8; t_{est} - t_{sim} [ns] ",
1324  100,
1325  -1.,
1326  1.);
1327  }
1328 
1329  meBarrelPIDp_ = ibook.book1D("BarrelPIDp", "PID track MTD momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
1330  meEndcapPIDp_ = ibook.book1D("EndcapPIDp", "PID track with MTD momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
1331 
1332  meBarrelNoPIDtype_ = ibook.book1D("BarrelNoPIDtype", "Barrel PID failure category", 4, 0., 4.);
1333  meEndcapNoPIDtype_ = ibook.book1D("EndcapNoPIDtype", "Endcap PID failure category", 4, 0., 4.);
1334 
1335  meBarrelTruePiNoPID_ =
1336  ibook.book1D("BarrelTruePiNoPID", "True pi NoPID momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
1337  meBarrelTrueKNoPID_ =
1338  ibook.book1D("BarrelTrueKNoPID", "True k NoPID momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
1339  meBarrelTruePNoPID_ =
1340  ibook.book1D("BarrelTruePNoPID", "True p NoPID momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
1341  meEndcapTruePiNoPID_ =
1342  ibook.book1D("EndcapTruePiNoPID", "True pi NoPID momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
1343  meEndcapTrueKNoPID_ =
1344  ibook.book1D("EndcapTrueKNoPID", "True k NoPID momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
1345  meEndcapTruePNoPID_ =
1346  ibook.book1D("EndcapTruePNoPID", "True p NoPID momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
1347 
1348  meBarrelTruePiAsPi_ =
1349  ibook.book1D("BarrelTruePiAsPi", "True pi as pi momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
1350  meBarrelTruePiAsK_ =
1351  ibook.book1D("BarrelTruePiAsK", "True pi as k momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
1352  meBarrelTruePiAsP_ =
1353  ibook.book1D("BarrelTruePiAsP", "True pi as p momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
1354  meEndcapTruePiAsPi_ =
1355  ibook.book1D("EndcapTruePiAsPi", "True pi as pi momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
1356  meEndcapTruePiAsK_ =
1357  ibook.book1D("EndcapTruePiAsK", "True pi as k momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
1358  meEndcapTruePiAsP_ =
1359  ibook.book1D("EndcapTruePiAsP", "True pi as p momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
1360 
1361  meBarrelTrueKAsPi_ =
1362  ibook.book1D("BarrelTrueKAsPi", "True k as pi momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
1363  meBarrelTrueKAsK_ = ibook.book1D("BarrelTrueKAsK", "True k as k momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
1364  meBarrelTrueKAsP_ = ibook.book1D("BarrelTrueKAsP", "True k as p momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
1365  meEndcapTrueKAsPi_ =
1366  ibook.book1D("EndcapTrueKAsPi", "True k as pi momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
1367  meEndcapTrueKAsK_ = ibook.book1D("EndcapTrueKAsK", "True k as k momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
1368  meEndcapTrueKAsP_ = ibook.book1D("EndcapTrueKAsP", "True k as p momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
1369 
1370  meBarrelTruePAsPi_ =
1371  ibook.book1D("BarrelTruePAsPi", "True p as pi momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
1372  meBarrelTruePAsK_ = ibook.book1D("BarrelTruePAsK", "True p as k momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
1373  meBarrelTruePAsP_ = ibook.book1D("BarrelTruePAsP", "True p as p momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
1374  meEndcapTruePAsPi_ =
1375  ibook.book1D("EndcapTruePAsPi", "True p as pi momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
1376  meEndcapTruePAsK_ = ibook.book1D("EndcapTruePAsK", "True p as k momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
1377  meEndcapTruePAsP_ = ibook.book1D("EndcapTruePAsP", "True p as p momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
1378 }
1379 
1380 bool Primary4DVertexValidation::matchRecoTrack2SimSignal(const reco::TrackBaseRef& recoTrack) {
1381  auto found = r2s_->find(recoTrack);
1382 
1383  // reco track not matched to any TP
1384  if (found == r2s_->end())
1385  return false;
1386 
1388  for (const auto& tp : found->val) {
1389  if (tp.first->eventId().bunchCrossing() == 0 && tp.first->eventId().event() == 0)
1390  return true;
1391  }
1392 
1393  // reco track not matched to any TP from signal vertex
1394  return false;
1395 }
1396 
1397 std::pair<const edm::Ref<std::vector<TrackingParticle>>*, int> Primary4DVertexValidation::getMatchedTP(
1398  const reco::TrackBaseRef& recoTrack, const TrackingVertexRef& vsim) {
1399  auto found = r2s_->find(recoTrack);
1400 
1401  // reco track not matched to any TP (fake tracks)
1402  if (found == r2s_->end())
1403  return std::make_pair(nullptr, -1);
1404 
1405  // matched TP equal to any TP of a given sim vertex
1406  for (const auto& tp : found->val) {
1407  if (std::find_if(vsim->daughterTracks_begin(), vsim->daughterTracks_end(), [&](const TrackingParticleRef& vtp) {
1408  return tp.first == vtp;
1409  }) != vsim->daughterTracks_end())
1410  return std::make_pair(&tp.first, 0);
1411  // matched TP not associated to any daughter track of a given sim vertex but having the same eventID (track from secondary vtx)
1412  else if (tp.first->eventId().bunchCrossing() == vsim->eventId().bunchCrossing() &&
1413  tp.first->eventId().event() == vsim->eventId().event()) {
1414  return std::make_pair(&tp.first, 1);
1415  }
1416  // matched TP not associated to any sim vertex of a given simulated event (PU track)
1417  else {
1418  return std::make_pair(&tp.first, 2);
1419  }
1420  }
1421 
1422  // reco track not matched to any TP from vertex
1423  return std::make_pair(nullptr, -1);
1424 }
1425 
1426 double Primary4DVertexValidation::timeFromTrueMass(double mass, double pathlength, double momentum, double time) {
1427  if (time > 0 && pathlength > 0 && mass > 0) {
1428  double gammasq = 1. + momentum * momentum / (mass * mass);
1429  double v = c_ * std::sqrt(1. - 1. / gammasq); // cm / ns
1430  double t_est = time - (pathlength / v);
1431 
1432  return t_est;
1433  } else {
1434  return -1;
1435  }
1436 }
1437 
1438 bool Primary4DVertexValidation::select(const reco::Vertex& v, int level) {
1439  /* level
1440  0 !isFake && ndof>4 (default)
1441  1 !isFake && ndof>4 && prob > 0.01
1442  2 !isFake && ndof>4 && prob > 0.01 && ptmax2 > 0.4
1443  */
1444  if (v.isFake())
1445  return false;
1446  if ((level == 0) && (v.ndof() > selNdof_))
1447  return true;
1448  /*if ((level == 1) && (v.ndof() > selNdof_) && (vertex_pxy(v) > 0.01))
1449  return true;
1450  if ((level == 2) && (v.ndof() > selNdof_) && (vertex_pxy(v) > 0.01) && (vertex_ptmax2(v) > 0.4))
1451  return true;
1452  if ((level == 3) && (v.ndof() > selNdof_) && (vertex_ptmax2(v) < 0.4))
1453  return true;*/
1454  return false;
1455 }
1456 
1457 void Primary4DVertexValidation::observablesFromJets(const std::vector<reco::Track>& reco_Tracks,
1458  const std::vector<double>& mass_Tracks,
1459  const std::vector<int>& category_Tracks,
1460  const std::string& skip_Tracks,
1461  unsigned int& n_Jets,
1462  double& sum_EtJets,
1463  double& sum_Pt2Jets,
1464  double& met_Pt,
1465  double& sum_PzJets) {
1466  double sum_PtJets = 0;
1467  n_Jets = 0;
1468  sum_EtJets = 0;
1469  sum_Pt2Jets = 0;
1470  met_Pt = 0;
1471  sum_PzJets = 0;
1472  auto met = LorentzVector(0, 0, 0, 0);
1473  std::vector<fastjet::PseudoJet> fjInputs_;
1474  fjInputs_.clear();
1475  size_t countScale0 = 0;
1476  for (size_t i = 0; i < reco_Tracks.size(); i++) {
1477  const auto recotr = reco_Tracks[i];
1478  const auto mass = mass_Tracks[i];
1479  float scale = 1.;
1480  if (recotr.charge() == 0) {
1481  continue;
1482  }
1483  // skip PU tracks in jet definition if skip_PU is required
1484  if (skip_Tracks == "skip_PU" && category_Tracks[i] == 2) {
1485  continue;
1486  }
1487  // skip fake tracks in jet definition if skip_Fake is required
1488  if (skip_Tracks == "skip_Fake" && category_Tracks[i] == -1) {
1489  continue;
1490  }
1491  if (recotr.pt() != 0) {
1492  scale = (recotr.pt() - recotr.ptError()) / recotr.pt();
1493  }
1494  if (edm::isNotFinite(scale)) {
1495  edm::LogWarning("Primary4DVertexValidation") << "Scaling is NAN ignoring this recotrack" << std::endl;
1496  scale = 0;
1497  }
1498  if (scale < 0) {
1499  scale = 0;
1500  countScale0++;
1501  }
1502  if (scale != 0) {
1503  fjInputs_.push_back(fastjet::PseudoJet(recotr.px() * scale,
1504  recotr.py() * scale,
1505  recotr.pz() * scale,
1506  std::sqrt(recotr.p() * recotr.p() + mass * mass) * scale));
1507  }
1508  }
1509  fastjet::ClusterSequence sequence(fjInputs_, fastjet::JetDefinition(fastjet::antikt_algorithm, 0.4));
1510  auto jets = fastjet::sorted_by_pt(sequence.inclusive_jets(0));
1511  for (const auto& pj : jets) {
1512  auto p4 = LorentzVector(pj.px(), pj.py(), pj.pz(), pj.e());
1513  sum_EtJets += std::sqrt(p4.e() * p4.e() - p4.P() * p4.P() + p4.pt() * p4.pt());
1514  sum_PtJets += p4.pt();
1515  sum_Pt2Jets += (p4.pt() * p4.pt() * 0.8 * 0.8);
1516  met += p4;
1517  sum_PzJets += p4.pz();
1518  n_Jets++;
1519  }
1520  met_Pt = met.pt();
1521  double metAbove = met_Pt - 2 * std::sqrt(sum_PtJets);
1522  if (metAbove > 0) {
1523  sum_Pt2Jets += (metAbove * metAbove);
1524  }
1525  if (countScale0 == reco_Tracks.size()) {
1526  sum_Pt2Jets = countScale0 * 0.01; //leave some epsilon value to sort vertices with unknown pt
1527  }
1528 }
1529 
1530 void Primary4DVertexValidation::isParticle(const reco::TrackBaseRef& recoTrack,
1531  const edm::ValueMap<float>& sigmat0,
1532  const edm::ValueMap<float>& sigmat0Safe,
1533  const edm::ValueMap<float>& probPi,
1534  const edm::ValueMap<float>& probK,
1535  const edm::ValueMap<float>& probP,
1536  unsigned int& no_PIDtype,
1537  bool& no_PID,
1538  bool& is_Pi,
1539  bool& is_K,
1540  bool& is_P) {
1541  no_PIDtype = 0;
1542  no_PID = false;
1543  is_Pi = false;
1544  is_K = false;
1545  is_P = false;
1546  if (probPi[recoTrack] == -1) {
1547  no_PIDtype = 1;
1548  } else if (edm::isNotFinite(probPi[recoTrack])) {
1549  no_PIDtype = 2;
1550  } else if (probPi[recoTrack] == 1 && probK[recoTrack] == 0 && probP[recoTrack] == 0 &&
1551  sigmat0[recoTrack] < sigmat0Safe[recoTrack]) {
1552  no_PIDtype = 3;
1553  }
1554  no_PID = no_PIDtype > 0;
1555  is_Pi = !no_PID && 1. - probPi[recoTrack] < minProbHeavy_;
1556  is_K = !no_PID && !is_Pi && probK[recoTrack] > probP[recoTrack];
1557  is_P = !no_PID && !is_Pi && !is_K;
1558 }
1559 
1560 void Primary4DVertexValidation::getWosWnt(const reco::Vertex& recoVtx,
1561  const reco::TrackBaseRef& recoTrk,
1562  const edm::ValueMap<float>& sigmat0,
1563  const edm::ValueMap<float>& mtdQualMVA,
1564  const edm::Handle<reco::BeamSpot>& BS,
1565  double& wos,
1566  double& wnt) {
1567  double dz2_beam = pow((*BS).BeamWidthX() * cos(recoTrk->phi()) / tan(recoTrk->theta()), 2) +
1568  pow((*BS).BeamWidthY() * sin(recoTrk->phi()) / tan(recoTrk->theta()), 2);
1569  double dz2 =
1570  pow(recoTrk->dzError(), 2) + dz2_beam + pow(0.0020, 2); // added 20 um, some tracks have crazy small resolutions
1571  wos = recoVtx.trackWeight(recoTrk) / dz2;
1572  wnt = recoVtx.trackWeight(recoTrk) *
1573  std::min(recoTrk->pt(), 1.0); // pt-weighted number of tracks (downweights pt < 1 GeV tracks)
1574 
1575  // If tracks have time information, give more weight to tracks with good resolution
1576  if (sigmat0[recoTrk] > 0 && mtdQualMVA[recoTrk] > mvaTh_) {
1577  double sigmaZ = (*BS).sigmaZ();
1578  double sigmaT = sigmaZ / c_; // c in cm/ns
1579  wos = wos / erf(sigmat0[recoTrk] / sigmaT);
1580  }
1581 }
1582 
1583 /* Extract information form TrackingParticles/TrackingVertex and fill
1584  * the helper class simPrimaryVertex with proper generation-level
1585  * information */
1586 std::vector<Primary4DVertexValidation::simPrimaryVertex> Primary4DVertexValidation::getSimPVs(
1587  const edm::Handle<TrackingVertexCollection>& tVC) {
1588  std::vector<Primary4DVertexValidation::simPrimaryVertex> simpv;
1589  int current_event = -1;
1590  int s = -1;
1591  for (TrackingVertexCollection::const_iterator v = tVC->begin(); v != tVC->end(); ++v) {
1592  // We keep only the first vertex from all the events at BX=0.
1593  if (v->eventId().bunchCrossing() != 0)
1594  continue;
1595  if (v->eventId().event() != current_event) {
1596  current_event = v->eventId().event();
1597  } else {
1598  continue;
1599  }
1600  s++;
1601  if (std::abs(v->position().z()) > 1000)
1602  continue; // skip junk vertices
1603 
1604  // could be a new vertex, check all primaries found so far to avoid multiple entries
1605  simPrimaryVertex sv(v->position().x(), v->position().y(), v->position().z(), v->position().t());
1606  sv.eventId = v->eventId();
1607  sv.sim_vertex = TrackingVertexRef(tVC, std::distance(tVC->begin(), v));
1608  sv.OriginalIndex = s;
1609 
1610  for (TrackingParticleRefVector::iterator iTrack = v->daughterTracks_begin(); iTrack != v->daughterTracks_end();
1611  ++iTrack) {
1612  assert((**iTrack).eventId().bunchCrossing() == 0);
1613  }
1614  simPrimaryVertex* vp = nullptr; // will become non-NULL if a vertex is found and then point to it
1615  for (std::vector<simPrimaryVertex>::iterator v0 = simpv.begin(); v0 != simpv.end(); v0++) {
1616  if ((sv.eventId == v0->eventId) && (std::abs(sv.x - v0->x) < 1e-5) && (std::abs(sv.y - v0->y) < 1e-5) &&
1617  (std::abs(sv.z - v0->z) < 1e-5)) {
1618  vp = &(*v0);
1619  break;
1620  }
1621  }
1622  if (!vp) {
1623  // this is a new vertex, add it to the list of sim-vertices
1624  simpv.push_back(sv);
1625  vp = &simpv.back();
1626  }
1627 
1628  // Loop over daughter track(s) as Tracking Particles
1629  for (TrackingVertex::tp_iterator iTP = v->daughterTracks_begin(); iTP != v->daughterTracks_end(); ++iTP) {
1630  auto momentum = (*(*iTP)).momentum();
1631  const reco::Track* matched_best_reco_track = nullptr;
1632  double match_quality = -1;
1633  if (use_only_charged_tracks_ && (**iTP).charge() == 0)
1634  continue;
1635  if (s2r_->find(*iTP) != s2r_->end()) {
1636  matched_best_reco_track = (*s2r_)[*iTP][0].first.get();
1637  match_quality = (*s2r_)[*iTP][0].second;
1638  }
1639 
1640  vp->ptot.setPx(vp->ptot.x() + momentum.x());
1641  vp->ptot.setPy(vp->ptot.y() + momentum.y());
1642  vp->ptot.setPz(vp->ptot.z() + momentum.z());
1643  vp->ptot.setE(vp->ptot.e() + (**iTP).energy());
1644  vp->pt += (**iTP).pt();
1645  vp->ptsq += ((**iTP).pt() * (**iTP).pt());
1646  vp->nGenTrk++;
1647 
1648  if (matched_best_reco_track) {
1649  vp->num_matched_reco_tracks++;
1650  vp->average_match_quality += match_quality;
1651  }
1652  } // End of for loop on daughters sim-particles
1653  if (vp->num_matched_reco_tracks) {
1654  vp->average_match_quality /= static_cast<float>(vp->num_matched_reco_tracks);
1655  }
1656  LogTrace("Primary4DVertexValidation")
1657  << "average number of associated tracks: " << vp->num_matched_reco_tracks / static_cast<float>(vp->nGenTrk)
1658  << " with average quality: " << vp->average_match_quality;
1659  } // End of for loop on tracking vertices
1660 
1661  // In case of no simulated vertices, break here
1662  if (simpv.empty())
1663  return simpv;
1664 
1665  // Now compute the closest distance in z between all simulated vertex
1666  // first initialize
1667  auto prev_z = simpv.back().z;
1668  for (simPrimaryVertex& vsim : simpv) {
1669  vsim.closest_vertex_distance_z = std::abs(vsim.z - prev_z);
1670  prev_z = vsim.z;
1671  }
1672  // then calculate
1673  for (std::vector<simPrimaryVertex>::iterator vsim = simpv.begin(); vsim != simpv.end(); vsim++) {
1674  std::vector<simPrimaryVertex>::iterator vsim2 = vsim;
1675  vsim2++;
1676  for (; vsim2 != simpv.end(); vsim2++) {
1677  double distance = std::abs(vsim->z - vsim2->z);
1678  // need both to be complete
1679  vsim->closest_vertex_distance_z = std::min(vsim->closest_vertex_distance_z, distance);
1680  vsim2->closest_vertex_distance_z = std::min(vsim2->closest_vertex_distance_z, distance);
1681  }
1682  }
1683  return simpv;
1684 }
1685 
1686 /* Extract information form reco Vertex and fill the helper class
1687  * recoPrimaryVertex with proper reco-level information */
1688 std::vector<Primary4DVertexValidation::recoPrimaryVertex> Primary4DVertexValidation::getRecoPVs(
1689  const edm::Handle<edm::View<reco::Vertex>>& tVC) {
1690  std::vector<Primary4DVertexValidation::recoPrimaryVertex> recopv;
1691  int r = -1;
1692  for (auto v = tVC->begin(); v != tVC->end(); ++v) {
1693  r++;
1694  // Skip junk vertices
1695  if (std::abs(v->z()) > 1000)
1696  continue;
1697  if (v->isFake() || !v->isValid())
1698  continue;
1699 
1700  recoPrimaryVertex sv(v->position().x(), v->position().y(), v->position().z());
1701  sv.recVtx = &(*v);
1702  sv.recVtxRef = reco::VertexBaseRef(tVC, std::distance(tVC->begin(), v));
1703 
1704  sv.OriginalIndex = r;
1705  sv.ndof = v->ndof();
1706  // this is a new vertex, add it to the list of reco-vertices
1707  recopv.push_back(sv);
1708  Primary4DVertexValidation::recoPrimaryVertex* vp = &recopv.back();
1709 
1710  // Loop over daughter track(s)
1711  for (auto iTrack = v->tracks_begin(); iTrack != v->tracks_end(); ++iTrack) {
1712  auto momentum = (*(*iTrack)).innerMomentum();
1713  if (momentum.mag2() == 0)
1714  momentum = (*(*iTrack)).momentum();
1715  vp->pt += std::sqrt(momentum.perp2());
1716  vp->ptsq += (momentum.perp2());
1717  vp->nRecoTrk++;
1718 
1719  auto matched = r2s_->find(*iTrack);
1720  if (matched != r2s_->end()) {
1721  vp->num_matched_sim_tracks++;
1722  }
1723 
1724  } // End of for loop on daughters reconstructed tracks
1725  } // End of for loop on tracking vertices
1726 
1727  // In case of no reco vertices, break here
1728  if (recopv.empty())
1729  return recopv;
1730 
1731  // Now compute the closest distance in z between all reconstructed vertex
1732  // first initialize
1733  auto prev_z = recopv.back().z;
1734  for (recoPrimaryVertex& vreco : recopv) {
1735  vreco.closest_vertex_distance_z = std::abs(vreco.z - prev_z);
1736  prev_z = vreco.z;
1737  }
1738  for (std::vector<recoPrimaryVertex>::iterator vreco = recopv.begin(); vreco != recopv.end(); vreco++) {
1739  std::vector<recoPrimaryVertex>::iterator vreco2 = vreco;
1740  vreco2++;
1741  for (; vreco2 != recopv.end(); vreco2++) {
1742  double distance = std::abs(vreco->z - vreco2->z);
1743  // need both to be complete
1744  vreco->closest_vertex_distance_z = std::min(vreco->closest_vertex_distance_z, distance);
1745  vreco2->closest_vertex_distance_z = std::min(vreco2->closest_vertex_distance_z, distance);
1746  }
1747  }
1748  return recopv;
1749 }
1750 
1751 // ------------ method called to produce the data ------------
1752 void Primary4DVertexValidation::matchReco2Sim(std::vector<recoPrimaryVertex>& recopv,
1753  std::vector<simPrimaryVertex>& simpv,
1754  const edm::ValueMap<float>& sigmat0,
1755  const edm::ValueMap<float>& MVA,
1756  const edm::Handle<reco::BeamSpot>& BS) {
1757  // Initialization (clear wos and wnt)
1758  for (auto vv : simpv) {
1759  vv.wnt.clear();
1760  vv.wos.clear();
1761  }
1762  for (auto rv : recopv) {
1763  rv.wnt.clear();
1764  rv.wos.clear();
1765  }
1766 
1767  // Filling infos for matching rec and sim vertices
1768  for (unsigned int iv = 0; iv < recopv.size(); iv++) {
1769  const reco::Vertex* vertex = recopv.at(iv).recVtx;
1770  LogTrace("Primary4DVertexValidation") << "iv (rec): " << iv;
1771 
1772  for (unsigned int iev = 0; iev < simpv.size(); iev++) {
1773  double wnt = 0;
1774  double wos = 0;
1775  double evwnt = 0; // weighted number of tracks from sim event iev in the current recvtx
1776  double evwos = 0; // weight over sigma**2 of sim event iev in the current recvtx
1777  unsigned int evnt = 0; // number of tracks from sim event iev in the current recvtx
1778 
1779  for (auto iTrack = vertex->tracks_begin(); iTrack != vertex->tracks_end(); ++iTrack) {
1780  if (vertex->trackWeight(*iTrack) < trackweightTh_) {
1781  continue;
1782  }
1783 
1784  auto tp_info = getMatchedTP(*iTrack, simpv.at(iev).sim_vertex).first;
1785  int matchCategory = getMatchedTP(*iTrack, simpv.at(iev).sim_vertex).second;
1786  // matched TP equal to any TP of a given sim vertex
1787  if (tp_info != nullptr && matchCategory == 0) {
1788  getWosWnt(*vertex, *iTrack, MVA, sigmat0, BS, wos, wnt);
1789  simpv.at(iev).addTrack(iv, wos, wnt);
1790  recopv.at(iv).addTrack(iev, wos, wnt);
1791  evwos += wos;
1792  evwnt += wnt;
1793  evnt++;
1794  }
1795  } // RecoTracks loop
1796 
1797  // require 2 tracks for a wos-match
1798  if ((evwos > 0) && (evwos > recopv.at(iv).maxwos) && (evnt > 1)) {
1799  recopv.at(iv).wosmatch = iev;
1800  recopv.at(iv).maxwos = evwos;
1801  recopv.at(iv).maxwosnt = evnt;
1802  LogTrace("Primary4DVertexValidation") << "dominating sim event (iev): " << iev << " evwos = " << evwos;
1803  }
1804 
1805  // weighted track counting match, require at least one track
1806  if ((evnt > 0) && (evwnt > recopv.at(iv).maxwnt)) {
1807  recopv.at(iv).wntmatch = iev;
1808  recopv.at(iv).maxwnt = evwnt;
1809  }
1810  } // TrackingVertex loop
1811  if (recopv.at(iv).maxwos > 0) {
1812  simpv.at(recopv.at(iv).wosmatch).wos_dominated_recv.push_back(iv);
1813  simpv.at(recopv.at(iv).wosmatch).nwosmatch++; // count the rec vertices dominated by a sim vertex using wos
1814  assert(iv < recopv.size());
1815  }
1816  LogTrace("Primary4DVertexValidation") << "largest contribution to wos: wosmatch (iev) = " << recopv.at(iv).wosmatch
1817  << " maxwos = " << recopv.at(iv).maxwos;
1818  if (recopv.at(iv).maxwnt > 0) {
1819  simpv.at(recopv.at(iv).wntmatch).nwntmatch++; // count the rec vertices dominated by a sim vertex using wnt
1820  }
1821  } // RecoPrimaryVertex
1822 
1823  // reset
1824  for (auto& vrec : recopv) {
1825  vrec.sim = NOT_MATCHED;
1826  vrec.matchQuality = 0;
1827  }
1828  unsigned int iev = 0;
1829  for (auto& vv : simpv) {
1830  LogTrace("Primary4DVertexValidation") << "iev (sim): " << iev;
1831  LogTrace("Primary4DVertexValidation") << "wos_dominated_recv.size: " << vv.wos_dominated_recv.size();
1832  for (unsigned int i = 0; i < vv.wos_dominated_recv.size(); i++) {
1833  auto recov = vv.wos_dominated_recv.at(i);
1834  LogTrace("Primary4DVertexValidation")
1835  << "index of reco vertex: " << recov << " that has a wos: " << vv.wos.at(recov) << " at position " << i;
1836  }
1837  vv.rec = NOT_MATCHED;
1838  vv.matchQuality = 0;
1839  iev++;
1840  }
1841  // after filling infos, goes for the sim-reco match
1842  // this tries a one-to-one match, taking simPV with highest wos if there are > 1 simPV candidates
1843  for (unsigned int rank = 1; rank < maxRank_; rank++) {
1844  for (unsigned int iev = 0; iev < simpv.size(); iev++) { //loop on SimPV
1845  if (simpv.at(iev).rec != NOT_MATCHED) {
1846  continue; // only sim vertices not already matched
1847  }
1848  if (simpv.at(iev).nwosmatch == 0) {
1849  continue; // the sim vertex does not dominate any reco vertex
1850  }
1851  if (simpv.at(iev).nwosmatch > rank) {
1852  continue; // start with sim vertices dominating one rec vertex (rank 1), then go with the ones dominating more
1853  }
1854  unsigned int iv = NOT_MATCHED; // select a rec vertex index
1855  for (unsigned int k = 0; k < simpv.at(iev).wos_dominated_recv.size(); k++) {
1856  unsigned int rec = simpv.at(iev).wos_dominated_recv.at(k); //candidate rec vertex index
1857  auto vrec = recopv.at(rec);
1858  if (vrec.sim != NOT_MATCHED) {
1859  continue; // already matched
1860  }
1861  if (std::abs(simpv.at(iev).z - vrec.z) > zWosMatchMax_) {
1862  continue; // insanely far away
1863  }
1864  if ((iv == NOT_MATCHED) || simpv.at(iev).wos.at(rec) > simpv.at(iev).wos.at(iv)) {
1865  iv = rec;
1866  }
1867  }
1868  // if a viable candidate was found, make the link
1869  if (iv !=
1870  NOT_MATCHED) { // if the rec vertex has already been associated is possible that iv remains NOT_MATCHED at this point
1871  recopv.at(iv).sim = iev;
1872  simpv.at(iev).rec = iv;
1873  recopv.at(iv).matchQuality = rank;
1874  simpv.at(iev).matchQuality = rank;
1875  }
1876  } // iev
1877  } // rank
1878 
1879  // Reco vertices that are not necessarily dominated by a sim vertex, or whose dominating sim-vertex
1880  // has been matched already to another overlapping reco vertex, can still be associated to a specific
1881  // sim vertex (without being classified as dominating).
1882  // In terms of fitting 1d-distributions, this corresponds to a small peak close to a much larger nearby peak
1883  unsigned int ntry = 0;
1884  while (ntry++ < maxTry_) {
1885  unsigned nmatch = 0;
1886  for (unsigned int iev = 0; iev < simpv.size(); iev++) {
1887  if ((simpv.at(iev).rec != NOT_MATCHED) || (simpv.at(iev).wos.empty())) {
1888  continue;
1889  }
1890  // find a rec vertex for the NOT_MATCHED sim vertex
1891  unsigned int rec = NOT_MATCHED;
1892  for (auto rv : simpv.at(iev).wos) {
1893  if ((rec == NOT_MATCHED) || (rv.second > simpv.at(iev).wos.at(rec))) {
1894  rec = rv.first;
1895  }
1896  }
1897 
1898  if (rec == NOT_MATCHED) { // try with wnt match
1899  for (auto rv : simpv.at(iev).wnt) {
1900  if ((rec == NOT_MATCHED) || (rv.second > simpv.at(iev).wnt.at(rec))) {
1901  rec = rv.first;
1902  }
1903  }
1904  }
1905 
1906  if (rec == NOT_MATCHED) {
1907  continue; // should not happen
1908  }
1909  if (recopv.at(rec).sim != NOT_MATCHED) {
1910  continue; // already gone
1911  }
1912 
1913  // check if the rec vertex can be matched
1914  unsigned int rec2sim = NOT_MATCHED;
1915  for (auto sv : recopv.at(rec).wos) {
1916  if (simpv.at(sv.first).rec != NOT_MATCHED) {
1917  continue; // already used
1918  }
1919  if ((rec2sim == NOT_MATCHED) || (sv.second > recopv.at(rec).wos.at(rec2sim))) {
1920  rec2sim = sv.first;
1921  }
1922  }
1923  if (iev == rec2sim) {
1924  // do the match and assign lowest quality (i.e. max rank)
1925  recopv.at(rec).sim = iev;
1926  recopv.at(rec).matchQuality = maxRank_;
1927  simpv.at(iev).rec = rec;
1928  simpv.at(iev).matchQuality = maxRank_;
1929  nmatch++;
1930  }
1931  } // sim loop
1932  if (nmatch == 0) {
1933  break;
1934  }
1935  } // ntry
1936 
1937 // Debugging
1938 #ifdef EDM_ML_DEBUG
1939  unsigned int nmatch_tot = 0, n_dzgtsz = 0;
1940  unsigned int n_rank1 = 0, n_rank2 = 0, n_rank3 = 0, n_rank8 = 0;
1941 
1942  for (unsigned int iev = 0; iev < simpv.size(); iev++) { //loop on SimPV
1943  if (simpv.at(iev).rec != NOT_MATCHED) {
1944  unsigned int rec = simpv.at(iev).rec;
1945  unsigned int wosmatch = recopv.at(rec).wosmatch;
1946  LogTrace("Primary4DVertexValidation")
1947  << "Final match: iev (sim) = " << std::setw(4) << iev << " sim.rec = " << std::setw(4) << rec
1948  << " rec.wosmatch = " << std::setw(5) << wosmatch << " dZ/sigmaZ = " << std::setw(6) << std::setprecision(2)
1949  << std::abs((recopv.at(rec).z - simpv.at(iev).z) / recopv.at(rec).recVtx->zError())
1950  << " match qual = " << std::setw(1) << recopv.at(rec).matchQuality;
1951  nmatch_tot++;
1952  if (std::abs((recopv.at(rec).z - simpv.at(iev).z) / recopv.at(rec).recVtx->zError()) > 3.) {
1953  n_dzgtsz++;
1954  }
1955  if (recopv.at(rec).matchQuality == 1) {
1956  n_rank1++;
1957  }
1958  if (recopv.at(rec).matchQuality == 2) {
1959  n_rank2++;
1960  }
1961  if (recopv.at(rec).matchQuality == 3) {
1962  n_rank3++;
1963  }
1964  if (recopv.at(rec).matchQuality == 8) {
1965  n_rank8++;
1966  }
1967  }
1968  }
1969  LogTrace("Primary4DVertexValidation") << "n_sim = " << simpv.size() << " n_rec = " << recopv.size()
1970  << " nmatch_tot = " << nmatch_tot << " n(dZ>sigmaZ) = " << n_dzgtsz
1971  << " n_rank1 = " << n_rank1 << " n_rank2 = " << n_rank2
1972  << " n_rank3 = " << n_rank3 << " n_rank8 = " << n_rank8;
1973 #endif
1974 }
1975 
1976 void Primary4DVertexValidation::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
1977  using edm::Handle;
1978  using edm::View;
1979  using std::cout;
1980  using std::endl;
1981  using std::vector;
1982  using namespace reco;
1983 
1984  std::vector<float> pileUpInfo_z;
1985 
1986  // get the pileup information
1987  edm::Handle<std::vector<PileupSummaryInfo>> puinfoH;
1988  if (iEvent.getByToken(vecPileupSummaryInfoToken_, puinfoH)) {
1989  for (auto const& pu_info : *puinfoH.product()) {
1990  if (pu_info.getBunchCrossing() == 0) {
1991  pileUpInfo_z = pu_info.getPU_zpositions();
1992  break;
1993  }
1994  }
1995  }
1996 
1997  edm::Handle<TrackingParticleCollection> TPCollectionH;
1998  iEvent.getByToken(trackingParticleCollectionToken_, TPCollectionH);
1999  if (!TPCollectionH.isValid())
2000  edm::LogWarning("Primary4DVertexValidation") << "TPCollectionH is not valid";
2001 
2002  edm::Handle<TrackingVertexCollection> TVCollectionH;
2003  iEvent.getByToken(trackingVertexCollectionToken_, TVCollectionH);
2004  if (!TVCollectionH.isValid())
2005  edm::LogWarning("Primary4DVertexValidation") << "TVCollectionH is not valid";
2006 
2007  edm::Handle<reco::SimToRecoCollection> simToRecoH;
2008  iEvent.getByToken(simToRecoAssociationToken_, simToRecoH);
2009  if (simToRecoH.isValid())
2010  s2r_ = simToRecoH.product();
2011  else
2012  edm::LogWarning("Primary4DVertexValidation") << "simToRecoH is not valid";
2013 
2014  edm::Handle<reco::RecoToSimCollection> recoToSimH;
2015  iEvent.getByToken(recoToSimAssociationToken_, recoToSimH);
2016  if (recoToSimH.isValid())
2017  r2s_ = recoToSimH.product();
2018  else
2019  edm::LogWarning("Primary4DVertexValidation") << "recoToSimH is not valid";
2020 
2021  edm::Handle<reco::BeamSpot> BeamSpotH;
2022  iEvent.getByToken(RecBeamSpotToken_, BeamSpotH);
2023  if (!BeamSpotH.isValid())
2024  edm::LogWarning("Primary4DVertexValidation") << "BeamSpotH is not valid";
2025 
2026  std::vector<simPrimaryVertex> simpv; // a list of simulated primary MC vertices
2027  simpv = getSimPVs(TVCollectionH);
2028  // this bool check if first vertex in that with highest pT
2029  bool signal_is_highest_pt =
2030  std::max_element(simpv.begin(), simpv.end(), [](const simPrimaryVertex& lhs, const simPrimaryVertex& rhs) {
2031  return lhs.ptsq < rhs.ptsq;
2032  }) == simpv.begin();
2033 
2034  std::vector<recoPrimaryVertex> recopv; // a list of reconstructed primary MC vertices
2035  edm::Handle<edm::View<reco::Vertex>> recVtxs;
2036  iEvent.getByToken(Rec4DVerToken_, recVtxs);
2037  if (!recVtxs.isValid())
2038  edm::LogWarning("Primary4DVertexValidation") << "recVtxs is not valid";
2039  recopv = getRecoPVs(recVtxs);
2040 
2041  const auto& trackAssoc = iEvent.get(trackAssocToken_);
2042  const auto& pathLength = iEvent.get(pathLengthToken_);
2043  const auto& momentum = iEvent.get(momentumToken_);
2044  const auto& time = iEvent.get(timeToken_);
2045  const auto& t0Pid = iEvent.get(t0PidToken_);
2046  const auto& sigmat0 = iEvent.get(sigmat0PidToken_);
2047  const auto& t0Safe = iEvent.get(t0SafePidToken_);
2048  const auto& sigmat0Safe = iEvent.get(sigmat0SafePidToken_);
2049  const auto& mtdQualMVA = iEvent.get(trackMVAQualToken_);
2050  const auto& tMtd = iEvent.get(tmtdToken_);
2051  const auto& tofPi = iEvent.get(tofPiToken_);
2052  const auto& tofK = iEvent.get(tofKToken_);
2053  const auto& tofP = iEvent.get(tofPToken_);
2054  const auto& probPi = iEvent.get(probPiToken_);
2055  const auto& probK = iEvent.get(probKToken_);
2056  const auto& probP = iEvent.get(probPToken_);
2057  const auto& fPDGTable = iSetup.getHandle(pdtToken_);
2058 
2059  // I have simPV and recoPV collections
2060  matchReco2Sim(recopv, simpv, sigmat0Safe, mtdQualMVA, BeamSpotH);
2061 
2062  // Loop on tracks
2063  for (unsigned int iv = 0; iv < recopv.size(); iv++) {
2064  if (recopv.at(iv).ndof > selNdof_) {
2065  const reco::Vertex* vertex = recopv.at(iv).recVtx;
2066 
2067  for (unsigned int iev = 0; iev < simpv.size(); iev++) {
2068  auto vsim = simpv.at(iev).sim_vertex;
2069 
2070  bool selectedVtxMatching = recopv.at(iv).sim == iev && simpv.at(iev).rec == iv;
2071  bool selectedLV = simpv.at(iev).eventId.bunchCrossing() == 0 && simpv.at(iev).eventId.event() == 0 &&
2072  recopv.at(iv).OriginalIndex == 0;
2073  bool selectedLVMatching = selectedVtxMatching && selectedLV; // bool for reco vtx leading match
2074  if (selectedLVMatching && !recopv.at(iv).is_signal()) {
2075  edm::LogWarning("Primary4DVertexValidation")
2076  << "Reco vtx leading match inconsistent: BX/ID " << simpv.at(iev).eventId.bunchCrossing() << " "
2077  << simpv.at(iev).eventId.event();
2078  }
2079 #ifdef EDM_ML_DEBUG
2080  if (selectedLVMatching) {
2081  printSimVtxRecoVtxInfo(simpv.at(iev), recopv.at(iv));
2082  }
2083 #endif
2084  double vzsim = simpv.at(iev).z;
2085  double vtsim = simpv.at(iev).t * simUnit_;
2086 
2087  double wnt = 0, wos = 0;
2088  double PUsumWnt = 0, PUsumWos = 0, SecsumWos = 0, FakesumWos = 0, PUsumPt = 0, PUsumPt2 = 0;
2089  double sumWnt = 0, sumWos = 0, sumPt = 0, sumPt2 = 0;
2090  unsigned int nt = 0, PUnt = 0, Fakent = 0;
2091 
2092  std::vector<double> massVector;
2093  std::vector<reco::Track> recotracks;
2094  std::vector<int> categoryVector;
2095  double sumEtJets = 0, sumPt2Jets = 0, metPt = 0, sumPzJets = 0;
2096  double sumEtJetsnoPU = 0, sumPt2JetsnoPU = 0, metPtnoPU = 0, sumPzJetsnoPU = 0;
2097  double sumEtJetsnoFake = 0, sumPt2JetsnoFake = 0, metPtnoFake = 0, sumPzJetsnoFake = 0;
2098  unsigned int nJets = 0, nJetsnoPU = 0, nJetsnoFake = 0;
2099  for (auto iTrack = vertex->tracks_begin(); iTrack != vertex->tracks_end(); ++iTrack) {
2100  if (trackAssoc[*iTrack] == -1) {
2101  edm::LogWarning("mtdTracks") << "Extended track not associated";
2102  continue;
2103  }
2104 
2105  // monitor all track weights associated to a vertex before selection on it
2106  if (selectedVtxMatching) {
2107  meVtxTrackW_->Fill(vertex->trackWeight(*iTrack));
2108  if (selectedLV) {
2109  meVtxTrackRecLVW_->Fill(vertex->trackWeight(*iTrack));
2110  }
2111  }
2112 
2113  if (vertex->trackWeight(*iTrack) < trackweightTh_)
2114  continue;
2115  bool noCrack = std::abs((*iTrack)->eta()) < trackMaxBtlEta_ || std::abs((*iTrack)->eta()) > trackMinEtlEta_;
2116 
2117  bool selectRecoTrk = trkRecSel(**iTrack);
2118  if (selectedLVMatching && selectRecoTrk) {
2119  if (noCrack) {
2120  meTrackEffPtTot_->Fill((*iTrack)->pt());
2121  }
2122  meTrackEffEtaTot_->Fill(std::abs((*iTrack)->eta()));
2123  }
2124 
2125  auto tp_info = getMatchedTP(*iTrack, vsim).first;
2126  int matchCategory = getMatchedTP(*iTrack, vsim).second;
2127 
2128  // PU, fake and secondary tracks
2129  if (selectedVtxMatching) {
2130  unsigned int no_PIDtype = 0;
2131  bool no_PID, is_Pi, is_K, is_P;
2132  int PartID = 211; // pion
2133  isParticle(*iTrack, sigmat0, sigmat0Safe, probPi, probK, probP, no_PIDtype, no_PID, is_Pi, is_K, is_P);
2134  if (!use3dNoTime_) {
2135  if (no_PID || is_Pi) {
2136  PartID = 211;
2137  } else if (is_K) {
2138  PartID = 321;
2139  } else if (is_P) {
2140  PartID = 2212;
2141  }
2142  }
2143  const HepPDT::ParticleData* PData = fPDGTable->particle(HepPDT::ParticleID(abs(PartID)));
2144  double mass = PData->mass().value();
2145  massVector.push_back(mass);
2146  recotracks.push_back(**iTrack);
2147  getWosWnt(*vertex, *iTrack, sigmat0Safe, mtdQualMVA, BeamSpotH, wos, wnt);
2148  meVtxTrackWnt_->Fill(wnt);
2149  if (selectedLV) {
2150  meVtxTrackRecLVWnt_->Fill(wnt);
2151  }
2152  // reco track matched to any TP
2153  if (tp_info != nullptr) {
2154 #ifdef EDM_ML_DEBUG
2155  if (selectedLV) {
2156  printMatchedRecoTrackInfo(*vertex, *iTrack, *tp_info, matchCategory);
2157  }
2158 #endif
2159  // matched TP not associated to any daughter track of a given sim vertex but having the same eventID (track from secondary vtx)
2160  if (matchCategory == 1) {
2161  categoryVector.push_back(matchCategory);
2162  SecsumWos += wos;
2163  }
2164  // matched TP not associated to any sim vertex of a given simulated event (PU track)
2165  if (matchCategory == 2) {
2166  if (optionalPlots_) {
2167  mePUTrackWnt_->Fill(wnt);
2168  if (selectedLV) {
2169  mePUTrackRecLVWnt_->Fill(wnt);
2170  }
2171  }
2172  PUsumWnt += wnt;
2173  PUsumWos += wos;
2174  PUsumPt += (*iTrack)->pt();
2175  PUsumPt2 += ((*iTrack)->pt() * (*iTrack)->pt());
2176  PUnt++;
2177  categoryVector.push_back(2);
2178  }
2179  }
2180  // reco track not matched to any TP (fake tracks)
2181  else {
2182  categoryVector.push_back(matchCategory);
2183  FakesumWos += wos;
2184  Fakent++;
2185  }
2186  nt++;
2187  sumWnt += wnt;
2188  sumWos += wos;
2189  sumPt += (*iTrack)->pt();
2190  sumPt2 += ((*iTrack)->pt() * (*iTrack)->pt());
2191  }
2192 
2193  // matched TP equal to any TP of a given sim vertex
2194  if (tp_info != nullptr && matchCategory == 0) {
2195  categoryVector.push_back(matchCategory);
2196  double mass = (*tp_info)->mass();
2197  double tsim = (*tp_info)->parentVertex()->position().t() * simUnit_;
2198  double tEst = timeFromTrueMass(mass, pathLength[*iTrack], momentum[*iTrack], time[*iTrack]);
2199 
2200  double xsim = (*tp_info)->parentVertex()->position().x();
2201  double ysim = (*tp_info)->parentVertex()->position().y();
2202  double zsim = (*tp_info)->parentVertex()->position().z();
2203  double xPCA = (*iTrack)->vx();
2204  double yPCA = (*iTrack)->vy();
2205  double zPCA = (*iTrack)->vz();
2206 
2207  double dZ = zPCA - zsim;
2208  double d3D = std::sqrt((xPCA - xsim) * (xPCA - xsim) + (yPCA - ysim) * (yPCA - ysim) + dZ * dZ);
2209  // orient d3D according to the projection of RECO - SIM onto simulated momentum
2210  if ((xPCA - xsim) * ((*tp_info)->px()) + (yPCA - ysim) * ((*tp_info)->py()) + dZ * ((*tp_info)->pz()) <
2211  0.) {
2212  d3D = -d3D;
2213  }
2214 
2215  // select TPs associated to the signal event
2216  bool selectTP = trkTPSelLV(**tp_info);
2217 
2218  if (selectedLVMatching && selectRecoTrk && selectTP) {
2219  meTrackMatchedTPZposResTot_->Fill((*iTrack)->vz() - vzsim);
2220  if (noCrack) {
2221  meTrackMatchedTPEffPtTot_->Fill((*iTrack)->pt());
2222  }
2223  meTrackMatchedTPEffEtaTot_->Fill(std::abs((*iTrack)->eta()));
2224  }
2225 
2226  if (sigmat0Safe[*iTrack] == -1)
2227  continue;
2228 
2229  if (selectedLVMatching && selectRecoTrk && selectTP) {
2230  meTrackMatchedTPResTot_->Fill(t0Safe[*iTrack] - vtsim);
2231  meTrackMatchedTPPullTot_->Fill((t0Safe[*iTrack] - vtsim) / sigmat0Safe[*iTrack]);
2232  if (noCrack) {
2233  meTrackMatchedTPEffPtMtd_->Fill((*iTrack)->pt());
2234  }
2235  meTrackMatchedTPEffEtaMtd_->Fill(std::abs((*iTrack)->eta()));
2236 
2237  unsigned int noPIDtype = 0;
2238  bool noPID = false, isPi = false, isK = false, isP = false;
2239  isParticle(*iTrack, sigmat0, sigmat0Safe, probPi, probK, probP, noPIDtype, noPID, isPi, isK, isP);
2240 
2241  if ((isPi && std::abs(tMtd[*iTrack] - tofPi[*iTrack] - t0Pid[*iTrack]) > tol_) ||
2242  (isK && std::abs(tMtd[*iTrack] - tofK[*iTrack] - t0Pid[*iTrack]) > tol_) ||
2243  (isP && std::abs(tMtd[*iTrack] - tofP[*iTrack] - t0Pid[*iTrack]) > tol_)) {
2244  edm::LogWarning("Primary4DVertexValidation")
2245  << "No match between mass hyp. and time: " << std::abs((*tp_info)->pdgId()) << " mass hyp pi/k/p "
2246  << isPi << " " << isK << " " << isP << " t0/t0safe " << t0Pid[*iTrack] << " " << t0Safe[*iTrack]
2247  << " tMtd - tof pi/K/p " << tMtd[*iTrack] - tofPi[*iTrack] << " " << tMtd[*iTrack] - tofK[*iTrack]
2248  << " " << tMtd[*iTrack] - tofP[*iTrack] << " Prob pi/K/p " << probPi[*iTrack] << " "
2249  << probK[*iTrack] << " " << probP[*iTrack];
2250  }
2251 
2252  if (std::abs((*iTrack)->eta()) < trackMaxBtlEta_) {
2253  meBarrelPIDp_->Fill((*iTrack)->p());
2254  meBarrelNoPIDtype_->Fill(noPIDtype + 0.5);
2255  if (std::abs((*tp_info)->pdgId()) == 211) {
2256  if (noPID) {
2257  meBarrelTruePiNoPID_->Fill((*iTrack)->p());
2258  } else if (isPi) {
2259  meBarrelTruePiAsPi_->Fill((*iTrack)->p());
2260  } else if (isK) {
2261  meBarrelTruePiAsK_->Fill((*iTrack)->p());
2262  } else if (isP) {
2263  meBarrelTruePiAsP_->Fill((*iTrack)->p());
2264  } else {
2265  edm::LogWarning("Primary4DVertexValidation")
2266  << "No PID class: " << std::abs((*tp_info)->pdgId()) << " t0/t0safe " << t0Pid[*iTrack] << " "
2267  << t0Safe[*iTrack] << " Prob pi/K/p " << probPi[*iTrack] << " " << probK[*iTrack] << " "
2268  << probP[*iTrack];
2269  }
2270  } else if (std::abs((*tp_info)->pdgId()) == 321) {
2271  if (noPID) {
2272  meBarrelTrueKNoPID_->Fill((*iTrack)->p());
2273  } else if (isPi) {
2274  meBarrelTrueKAsPi_->Fill((*iTrack)->p());
2275  } else if (isK) {
2276  meBarrelTrueKAsK_->Fill((*iTrack)->p());
2277  } else if (isP) {
2278  meBarrelTrueKAsP_->Fill((*iTrack)->p());
2279  } else {
2280  edm::LogWarning("Primary4DVertexValidation")
2281  << "No PID class: " << std::abs((*tp_info)->pdgId()) << " t0/t0safe " << t0Pid[*iTrack] << " "
2282  << t0Safe[*iTrack] << " Prob pi/K/p " << probPi[*iTrack] << " " << probK[*iTrack] << " "
2283  << probP[*iTrack];
2284  }
2285  } else if (std::abs((*tp_info)->pdgId()) == 2212) {
2286  if (noPID) {
2287  meBarrelTruePNoPID_->Fill((*iTrack)->p());
2288  } else if (isPi) {
2289  meBarrelTruePAsPi_->Fill((*iTrack)->p());
2290  } else if (isK) {
2291  meBarrelTruePAsK_->Fill((*iTrack)->p());
2292  } else if (isP) {
2293  meBarrelTruePAsP_->Fill((*iTrack)->p());
2294  } else {
2295  edm::LogWarning("Primary4DVertexValidation")
2296  << "No PID class: " << std::abs((*tp_info)->pdgId()) << " t0/t0safe " << t0Pid[*iTrack] << " "
2297  << t0Safe[*iTrack] << " Prob pi/K/p " << probPi[*iTrack] << " " << probK[*iTrack] << " "
2298  << probP[*iTrack];
2299  }
2300  }
2301  } else if (std::abs((*iTrack)->eta()) > trackMinEtlEta_ && std::abs((*iTrack)->eta()) < trackMaxEtlEta_) {
2302  meEndcapPIDp_->Fill((*iTrack)->p());
2303  meEndcapNoPIDtype_->Fill(noPIDtype + 0.5);
2304  if (std::abs((*tp_info)->pdgId()) == 211) {
2305  if (noPID) {
2306  meEndcapTruePiNoPID_->Fill((*iTrack)->p());
2307  } else if (isPi) {
2308  meEndcapTruePiAsPi_->Fill((*iTrack)->p());
2309  } else if (isK) {
2310  meEndcapTruePiAsK_->Fill((*iTrack)->p());
2311  } else if (isP) {
2312  meEndcapTruePiAsP_->Fill((*iTrack)->p());
2313  } else {
2314  edm::LogWarning("Primary4DVertexValidation")
2315  << "No PID class: " << std::abs((*tp_info)->pdgId()) << " t0/t0safe " << t0Pid[*iTrack] << " "
2316  << t0Safe[*iTrack] << " Prob pi/K/p " << probPi[*iTrack] << " " << probK[*iTrack] << " "
2317  << probP[*iTrack];
2318  }
2319  } else if (std::abs((*tp_info)->pdgId()) == 321) {
2320  if (noPID) {
2321  meEndcapTrueKNoPID_->Fill((*iTrack)->p());
2322  } else if (isPi) {
2323  meEndcapTrueKAsPi_->Fill((*iTrack)->p());
2324  } else if (isK) {
2325  meEndcapTrueKAsK_->Fill((*iTrack)->p());
2326  } else if (isP) {
2327  meEndcapTrueKAsP_->Fill((*iTrack)->p());
2328  } else {
2329  edm::LogWarning("Primary4DVertexValidation")
2330  << "No PID class: " << std::abs((*tp_info)->pdgId()) << " t0/t0safe " << t0Pid[*iTrack] << " "
2331  << t0Safe[*iTrack] << " Prob pi/K/p " << probPi[*iTrack] << " " << probK[*iTrack] << " "
2332  << probP[*iTrack];
2333  }
2334  } else if (std::abs((*tp_info)->pdgId()) == 2212) {
2335  if (noPID) {
2336  meEndcapTruePNoPID_->Fill((*iTrack)->p());
2337  } else if (isPi) {
2338  meEndcapTruePAsPi_->Fill((*iTrack)->p());
2339  } else if (isK) {
2340  meEndcapTruePAsK_->Fill((*iTrack)->p());
2341  } else if (isP) {
2342  meEndcapTruePAsP_->Fill((*iTrack)->p());
2343  } else {
2344  edm::LogWarning("Primary4DVertexValidation")
2345  << "No PID class: " << std::abs((*tp_info)->pdgId()) << " t0/t0safe " << t0Pid[*iTrack] << " "
2346  << t0Safe[*iTrack] << " Prob pi/K/p " << probPi[*iTrack] << " " << probK[*iTrack] << " "
2347  << probP[*iTrack];
2348  }
2349  }
2350  }
2351  }
2352  meTrackResTot_->Fill(t0Safe[*iTrack] - tsim);
2353  meTrackPullTot_->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
2354  meTrackZposResTot_->Fill(dZ);
2355  if ((*iTrack)->p() <= 2) {
2356  meTrackResLowPTot_->Fill(t0Safe[*iTrack] - tsim);
2357  meTrackPullLowPTot_->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
2358  } else {
2359  meTrackResHighPTot_->Fill(t0Safe[*iTrack] - tsim);
2360  meTrackPullHighPTot_->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
2361  }
2362 
2363  if (mtdQualMVA[(*iTrack)] < mvaL_) {
2364  meTrackZposRes_[0]->Fill(dZ);
2365  meTrack3DposRes_[0]->Fill(d3D);
2366  meTrackRes_[0]->Fill(t0Safe[*iTrack] - tsim);
2367  meTrackPull_[0]->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
2368 
2369  if (optionalPlots_) {
2370  meTrackResMass_[0]->Fill(t0Safe[*iTrack] - tEst);
2371  meTrackResMassTrue_[0]->Fill(tEst - tsim);
2372  }
2373 
2374  if ((*iTrack)->p() <= 2) {
2375  meTrackResLowP_[0]->Fill(t0Safe[*iTrack] - tsim);
2376  meTrackPullLowP_[0]->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
2377  } else if ((*iTrack)->p() > 2) {
2378  meTrackResHighP_[0]->Fill(t0Safe[*iTrack] - tsim);
2379  meTrackPullHighP_[0]->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
2380  }
2381 
2382  if (optionalPlots_) {
2383  if (std::abs((*tp_info)->pdgId()) == 2212) {
2384  meTrackResMassProtons_[0]->Fill(t0Safe[*iTrack] - tEst);
2385  meTrackResMassTrueProtons_[0]->Fill(tEst - tsim);
2386  } else if (std::abs((*tp_info)->pdgId()) == 211) {
2387  meTrackResMassPions_[0]->Fill(t0Safe[*iTrack] - tEst);
2388  meTrackResMassTruePions_[0]->Fill(tEst - tsim);
2389  }
2390  }
2391 
2392  } else if (mtdQualMVA[(*iTrack)] > mvaL_ && mtdQualMVA[(*iTrack)] < mvaH_) {
2393  meTrackZposRes_[1]->Fill(dZ);
2394  meTrack3DposRes_[1]->Fill(d3D);
2395  meTrackRes_[1]->Fill(t0Safe[*iTrack] - tsim);
2396  meTrackPull_[1]->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
2397 
2398  if (optionalPlots_) {
2399  meTrackResMass_[1]->Fill(t0Safe[*iTrack] - tEst);
2400  meTrackResMassTrue_[1]->Fill(tEst - tsim);
2401  }
2402 
2403  if ((*iTrack)->p() <= 2) {
2404  meTrackResLowP_[1]->Fill(t0Safe[*iTrack] - tsim);
2405  meTrackPullLowP_[1]->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
2406  } else if ((*iTrack)->p() > 2) {
2407  meTrackResHighP_[1]->Fill(t0Safe[*iTrack] - tsim);
2408  meTrackPullHighP_[1]->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
2409  }
2410 
2411  if (optionalPlots_) {
2412  if (std::abs((*tp_info)->pdgId()) == 2212) {
2413  meTrackResMassProtons_[1]->Fill(t0Safe[*iTrack] - tEst);
2414  meTrackResMassTrueProtons_[1]->Fill(tEst - tsim);
2415  } else if (std::abs((*tp_info)->pdgId()) == 211) {
2416  meTrackResMassPions_[1]->Fill(t0Safe[*iTrack] - tEst);
2417  meTrackResMassTruePions_[1]->Fill(tEst - tsim);
2418  }
2419  }
2420 
2421  } else if (mtdQualMVA[(*iTrack)] > mvaH_) {
2422  meTrackZposRes_[2]->Fill(dZ);
2423  meTrack3DposRes_[2]->Fill(d3D);
2424  meTrackRes_[2]->Fill(t0Safe[*iTrack] - tsim);
2425  meTrackPull_[2]->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
2426 
2427  if (optionalPlots_) {
2428  meTrackResMass_[2]->Fill(t0Safe[*iTrack] - tEst);
2429  meTrackResMassTrue_[2]->Fill(tEst - tsim);
2430  }
2431 
2432  if ((*iTrack)->p() <= 2) {
2433  meTrackResLowP_[2]->Fill(t0Safe[*iTrack] - tsim);
2434  meTrackPullLowP_[2]->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
2435  } else if ((*iTrack)->p() > 2) {
2436  meTrackResHighP_[2]->Fill(t0Safe[*iTrack] - tsim);
2437  meTrackPullHighP_[2]->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
2438  }
2439 
2440  if (optionalPlots_) {
2441  if (std::abs((*tp_info)->pdgId()) == 2212) {
2442  meTrackResMassProtons_[2]->Fill(t0Safe[*iTrack] - tEst);
2443  meTrackResMassTrueProtons_[2]->Fill(tEst - tsim);
2444  } else if (std::abs((*tp_info)->pdgId()) == 211) {
2445  meTrackResMassPions_[2]->Fill(t0Safe[*iTrack] - tEst);
2446  meTrackResMassTruePions_[2]->Fill(tEst - tsim);
2447  }
2448  }
2449  }
2450  } // if tp_info != nullptr && MatchCategory == 0
2451  } // loop on reco tracks
2452  if (selectedVtxMatching) {
2453  meVtxTrackMult_->Fill(log10(nt));
2454  mePUTrackRelMult_->Fill(static_cast<double>(PUnt) / nt);
2455  meFakeTrackRelMult_->Fill(static_cast<double>(Fakent) / nt);
2456  mePUTrackRelSumWnt_->Fill(PUsumWnt / sumWnt);
2457  mePUTrackRelSumWos_->Fill(PUsumWos / sumWos);
2458  meSecTrackRelSumWos_->Fill(SecsumWos / sumWos);
2459  meFakeTrackRelSumWos_->Fill(FakesumWos / sumWos);
2460  mePUTrackRelSumPt_->Fill(PUsumPt / sumPt);
2461  mePUTrackRelSumPt2_->Fill(PUsumPt2 / sumPt2);
2462 
2463  observablesFromJets(
2464  recotracks, massVector, categoryVector, "use_allTracks", nJets, sumEtJets, sumPt2Jets, metPt, sumPzJets);
2465  observablesFromJets(recotracks,
2466  massVector,
2467  categoryVector,
2468  "skip_Fake",
2469  nJetsnoFake,
2470  sumEtJetsnoFake,
2471  sumPt2JetsnoFake,
2472  metPtnoFake,
2473  sumPzJetsnoFake);
2474  observablesFromJets(recotracks,
2475  massVector,
2476  categoryVector,
2477  "skip_PU",
2478  nJetsnoPU,
2479  sumEtJetsnoPU,
2480  sumPt2JetsnoPU,
2481  metPtnoPU,
2482  sumPzJetsnoPU);
2483 
2484  meJetsPURelMult_->Fill(static_cast<double>(nJets - nJetsnoPU) / nJets);
2485  meJetsPURelHt_->Fill((sumEtJets - sumEtJetsnoPU) / sumEtJets);
2486  meJetsPURelSumPt2_->Fill((sumPt2Jets - sumPt2JetsnoPU) / sumPt2Jets);
2487  meJetsFakeRelSumPt2_->Fill((sumPt2Jets - sumPt2JetsnoFake) / sumPt2Jets);
2488  meJetsPURelMetPt_->Fill((metPt - metPtnoPU) / metPt);
2489  meJetsPURelSumPz_->Fill((sumPzJets - sumPzJetsnoPU) / sumPzJets);
2490 
2491  if (optionalPlots_) {
2492  mePUTrackMult_->Fill(PUnt);
2493  mePUTrackSumWnt_->Fill(log10(std::max(minThrSumWnt_, PUsumWnt)));
2494  mePUTrackSumWos_->Fill(log10(std::max(minThrSumWos_, PUsumWos)));
2495  meSecTrackSumWos_->Fill(log10(std::max(minThrSumWos_, SecsumWos)));
2496  mePUTrackSumPt_->Fill(log10(std::max(minThrSumPt_, PUsumPt)));
2497  mePUTrackSumPt2_->Fill(log10(std::max(minThrSumPt2_, PUsumPt2)));
2498 
2499  mePUTrackRelMultvsMult_->Fill(nt, static_cast<double>(PUnt) / nt);
2500  meFakeTrackRelMultvsMult_->Fill(nt, static_cast<double>(Fakent) / nt);
2501  mePUTrackRelSumWntvsSumWnt_->Fill(log10(std::max(minThrSumWnt_, sumWnt)), PUsumWnt / sumWnt);
2502  mePUTrackRelSumWosvsSumWos_->Fill(log10(std::max(minThrSumWos_, sumWos)), PUsumWos / sumWos);
2503  meSecTrackRelSumWosvsSumWos_->Fill(log10(std::max(minThrSumWos_, sumWos)), SecsumWos / sumWos);
2504  meFakeTrackRelSumWosvsSumWos_->Fill(log10(std::max(minThrSumWos_, sumWos)), FakesumWos / sumWos);
2505  mePUTrackRelSumPtvsSumPt_->Fill(log10(std::max(minThrSumPt_, sumPt)), PUsumPt / sumPt);
2506  mePUTrackRelSumPt2vsSumPt2_->Fill(log10(std::max(minThrSumPt2_, sumPt2)), PUsumPt2 / sumPt2);
2507 
2508  meJetsPUMult_->Fill(nJets - nJetsnoPU);
2509  meJetsPUHt_->Fill(log10(std::max(minThrSumPt_, sumEtJets - sumEtJetsnoPU)));
2510  meJetsPUSumPt2_->Fill(log10(std::max(minThrSumPt2_, sumPt2Jets - sumPt2JetsnoPU)));
2511  meJetsPUMetPt_->Fill(log10(std::max(minThrMetPt_, metPt - metPtnoPU)));
2512  meJetsPUSumPz_->Fill(log10(std::max(minThrSumPz_, std::abs(sumPzJets - sumPzJetsnoPU))));
2513 
2514  meJetsPURelMultvsMult_->Fill(nJets, static_cast<double>(nJets - nJetsnoPU) / nJets);
2515  meJetsPURelHtvsHt_->Fill(log10(std::max(minThrSumPt_, sumEtJets)), (sumEtJets - sumEtJetsnoPU) / sumEtJets);
2516  meJetsPURelSumPt2vsSumPt2_->Fill(log10(std::max(minThrSumPt2_, sumPt2Jets)),
2517  (sumPt2Jets - sumPt2JetsnoPU) / sumPt2Jets);
2518  meJetsFakeRelSumPt2vsSumPt2_->Fill(log10(std::max(minThrSumPt2_, sumPt2Jets)),
2519  (sumPt2Jets - sumPt2JetsnoFake) / sumPt2Jets);
2520  meJetsPURelMetPtvsMetPt_->Fill(log10(std::max(minThrMetPt_, metPt)), (metPt - metPtnoPU) / metPt);
2521  meJetsPURelSumPzvsSumPz_->Fill(log10(std::max(minThrSumPz_, std::abs(sumPzJets))),
2522  (sumPzJets - sumPzJetsnoPU) / sumPzJets);
2523  }
2524  if (selectedLV) {
2525  meVtxTrackRecLVMult_->Fill(log10(nt));
2526  mePUTrackRecLVRelMult_->Fill(static_cast<double>(PUnt) / nt);
2527  meFakeTrackRecLVRelMult_->Fill(static_cast<double>(Fakent) / nt);
2528  mePUTrackRecLVRelSumWnt_->Fill(PUsumWnt / sumWnt);
2529  mePUTrackRecLVRelSumWos_->Fill(PUsumWos / sumWos);
2530  meSecTrackRecLVRelSumWos_->Fill(SecsumWos / sumWos);
2531  meFakeTrackRecLVRelSumWos_->Fill(FakesumWos / sumWos);
2532  mePUTrackRecLVRelSumPt_->Fill(PUsumPt / sumPt);
2533  mePUTrackRecLVRelSumPt2_->Fill(PUsumPt2 / sumPt2);
2534 
2535  meJetsRecLVPURelMult_->Fill(static_cast<double>(nJets - nJetsnoPU) / nJets);
2536  meJetsRecLVPURelHt_->Fill((sumEtJets - sumEtJetsnoPU) / sumEtJets);
2537  meJetsRecLVPURelSumPt2_->Fill((sumPt2Jets - sumPt2JetsnoPU) / sumPt2Jets);
2538  meJetsRecLVFakeRelSumPt2_->Fill((sumPt2Jets - sumPt2JetsnoFake) / sumPt2Jets);
2539  meJetsRecLVPURelMetPt_->Fill((metPt - metPtnoPU) / metPt);
2540  meJetsRecLVPURelSumPz_->Fill((sumPzJets - sumPzJetsnoPU) / sumPzJets);
2541 
2542  LogTrace("Primary4DVertexValidation")
2543  << "#PUTrks = " << PUnt << " #Trks = " << nt << " PURelMult = " << std::setprecision(3)
2544  << static_cast<double>(PUnt) / nt;
2545  LogTrace("Primary4DVertexValidation")
2546  << "PUsumWnt = " << std::setprecision(3) << PUsumWnt << " sumWnt = " << std::setprecision(3) << sumWnt
2547  << " PURelsumWnt = " << std::setprecision(3) << PUsumWnt / sumWnt;
2548  LogTrace("Primary4DVertexValidation")
2549  << "PUsumWos = " << std::setprecision(3) << PUsumWos << " sumWos = " << std::setprecision(3) << sumWos
2550  << " PURelsumWos = " << std::setprecision(3) << PUsumWos / sumWos;
2551  LogTrace("Primary4DVertexValidation")
2552  << "PuSumPt = " << std::setprecision(3) << PUsumPt << " SumPt = " << std::setprecision(4) << sumPt
2553  << " PURelSumPt = " << std::setprecision(3) << PUsumPt / sumPt;
2554  LogTrace("Primary4DVertexValidation")
2555  << "PuSumPt2 = " << std::setprecision(3) << PUsumPt2 << " SumPt2 = " << std::setprecision(4) << sumPt2
2556  << " PURelSumPt2 = " << std::setprecision(3) << PUsumPt2 / sumPt2;
2557  if (optionalPlots_) {
2558  mePUTrackRecLVMult_->Fill(PUnt);
2559  mePUTrackRecLVSumWnt_->Fill(log10(std::max(minThrSumWnt_, PUsumWnt)));
2560  mePUTrackRecLVSumWos_->Fill(log10(std::max(minThrSumWos_, PUsumWos)));
2561  meSecTrackRecLVSumWos_->Fill(log10(std::max(minThrSumWos_, PUsumWos)));
2562  mePUTrackRecLVSumPt_->Fill(log10(std::max(minThrSumPt_, PUsumPt)));
2563  mePUTrackRecLVSumPt2_->Fill(log10(std::max(minThrSumPt2_, PUsumPt2)));
2564 
2565  mePUTrackRecLVRelMultvsMult_->Fill(nt, static_cast<double>(PUnt) / nt);
2566  meFakeTrackRecLVRelMultvsMult_->Fill(nt, static_cast<double>(Fakent) / nt);
2567  mePUTrackRecLVRelSumWntvsSumWnt_->Fill(log10(std::max(minThrSumWnt_, sumWnt)), PUsumWnt / sumWnt);
2568  mePUTrackRecLVRelSumWosvsSumWos_->Fill(log10(std::max(minThrSumWos_, sumWos)), PUsumWos / sumWos);
2569  meSecTrackRecLVRelSumWosvsSumWos_->Fill(log10(std::max(minThrSumWos_, sumWos)), SecsumWos / sumWos);
2570  meFakeTrackRecLVRelSumWosvsSumWos_->Fill(log10(std::max(minThrSumWos_, sumWos)), FakesumWos / sumWos);
2571  mePUTrackRecLVRelSumPtvsSumPt_->Fill(log10(std::max(minThrSumPt_, sumPt)), PUsumPt / sumPt);
2572  mePUTrackRecLVRelSumPt2vsSumPt2_->Fill(log10(std::max(minThrSumPt2_, sumPt2)), PUsumPt2 / sumPt2);
2573 
2574  meJetsRecLVPUMult_->Fill(nJets - nJetsnoPU);
2575  meJetsRecLVPUHt_->Fill(log10(std::max(minThrSumPt_, sumEtJets - sumEtJetsnoPU)));
2576  meJetsRecLVPUSumPt2_->Fill(log10(std::max(minThrSumPt2_, sumPt2Jets - sumPt2JetsnoPU)));
2577  meJetsRecLVPUMetPt_->Fill(log10(std::max(minThrMetPt_, metPt - metPtnoPU)));
2578  meJetsRecLVPUSumPz_->Fill(log10(std::max(minThrSumPz_, std::abs(sumPzJets - sumPzJetsnoPU))));
2579 
2580  meJetsRecLVPURelMultvsMult_->Fill(nJets, static_cast<double>(nJets - nJetsnoPU) / nJets);
2581  meJetsRecLVPURelHtvsHt_->Fill(log10(std::max(minThrSumPt_, sumEtJets)),
2582  (sumEtJets - sumEtJetsnoPU) / sumEtJets);
2583  meJetsRecLVPURelSumPt2vsSumPt2_->Fill(log10(std::max(minThrSumPt2_, sumPt2Jets)),
2584  (sumPt2Jets - sumPt2JetsnoPU) / sumPt2Jets);
2585  meJetsRecLVFakeRelSumPt2vsSumPt2_->Fill(log10(std::max(minThrSumPt2_, sumPt2Jets)),
2586  (sumPt2Jets - sumPt2JetsnoFake) / sumPt2Jets);
2587  meJetsRecLVPURelMetPtvsMetPt_->Fill(log10(std::max(minThrMetPt_, metPt)), (metPt - metPtnoPU) / metPt);
2588  meJetsRecLVPURelSumPzvsSumPz_->Fill(log10(std::max(minThrSumPz_, std::abs(sumPzJets))),
2589  (sumPzJets - sumPzJetsnoPU) / sumPzJets);
2590  }
2591  }
2592  }
2593  } // loop on simpv
2594  } // ndof
2595  } // loop on recopv
2596 
2597  int real = 0;
2598  int fake = 0;
2599  int other_fake = 0;
2600  int split = 0;
2601 
2602  auto puLineDensity = [&](double z) {
2603  // gaussian parameterization of line density vs z, z in cm, parameters in mm
2604  double argl = (z * 10. - lineDensityPar_[1]) / lineDensityPar_[2];
2605  return lineDensityPar_[0] * exp(-0.5 * argl * argl);
2606  };
2607 
2608  meRecVerNumber_->Fill(recopv.size());
2609  for (unsigned int ir = 0; ir < recopv.size(); ir++) {
2610  if (recopv.at(ir).ndof > selNdof_) {
2611  meRecoVtxVsLineDensity_->Fill(puLineDensity(recopv.at(ir).z));
2612  meRecPVZ_->Fill(recopv.at(ir).z, 1. / puLineDensity(recopv.at(ir).z));
2613  if (recopv.at(ir).recVtx->tError() > 0.) {
2614  meRecPVT_->Fill(recopv.at(ir).recVtx->t());
2615  }
2616  LogTrace("Primary4DVertexValidation") << "************* IR: " << ir;
2617  LogTrace("Primary4DVertexValidation")
2618  << "z: " << recopv.at(ir).z << " corresponding to line density: " << puLineDensity(recopv.at(ir).z);
2619  LogTrace("Primary4DVertexValidation") << "is_real: " << recopv.at(ir).is_real();
2620  LogTrace("Primary4DVertexValidation") << "is_fake: " << recopv.at(ir).is_fake();
2621  LogTrace("Primary4DVertexValidation") << "is_signal: " << recopv.at(ir).is_signal();
2622  LogTrace("Primary4DVertexValidation") << "split_from: " << recopv.at(ir).split_from();
2623  LogTrace("Primary4DVertexValidation") << "other fake: " << recopv.at(ir).other_fake();
2624  if (recopv.at(ir).is_real()) {
2625  real++;
2626  }
2627  if (recopv.at(ir).is_fake()) {
2628  fake++;
2629  }
2630  if (recopv.at(ir).other_fake()) {
2631  other_fake++;
2632  }
2633  if (recopv.at(ir).split_from() != -1) {
2634  split++;
2635  }
2636  } // ndof
2637  }
2638 
2639  LogTrace("Primary4DVertexValidation") << "is_real: " << real;
2640  LogTrace("Primary4DVertexValidation") << "is_fake: " << fake;
2641  LogTrace("Primary4DVertexValidation") << "split_from: " << split;
2642  LogTrace("Primary4DVertexValidation") << "other fake: " << other_fake;
2643  mePUvsRealV_->Fill(simpv.size(), real);
2644  mePUvsFakeV_->Fill(simpv.size(), fake);
2645  mePUvsOtherFakeV_->Fill(simpv.size(), other_fake);
2646  mePUvsSplitV_->Fill(simpv.size(), split);
2647 
2648  // fill vertices histograms here in a new loop
2649  for (unsigned int is = 0; is < simpv.size(); is++) {
2650  // protect against particle guns with very displaced vertices
2651  if (edm::isNotFinite(1. / puLineDensity(simpv.at(is).z))) {
2652  continue;
2653  }
2654  meSimPVZ_->Fill(simpv.at(is).z, 1. / puLineDensity(simpv.at(is).z));
2655  if (is == 0 && optionalPlots_) {
2656  meSimPosInSimOrigCollection_->Fill(simpv.at(is).OriginalIndex);
2657  }
2658 
2659  if (simpv.at(is).rec == NOT_MATCHED) {
2660  LogTrace("Primary4DVertexValidation") << "sim vertex: " << is << " is not matched with any reco";
2661  continue;
2662  }
2663 
2664  for (unsigned int ir = 0; ir < recopv.size(); ir++) {
2665  if (recopv.at(ir).ndof > selNdof_) {
2666  if (recopv.at(ir).sim == is && simpv.at(is).rec == ir) {
2667  meTimeRes_->Fill(recopv.at(ir).recVtx->t() - simpv.at(is).t * simUnit_);
2668  meTimePull_->Fill((recopv.at(ir).recVtx->t() - simpv.at(is).t * simUnit_) / recopv.at(ir).recVtx->tError());
2669  meMatchQual_->Fill(recopv.at(ir).matchQuality - 0.5);
2670  if (ir == 0) { // signal vertex plots
2671  meTimeSignalRes_->Fill(recopv.at(ir).recVtx->t() - simpv.at(is).t * simUnit_);
2672  meTimeSignalPull_->Fill((recopv.at(ir).recVtx->t() - simpv.at(is).t * simUnit_) /
2673  recopv.at(ir).recVtx->tError());
2674  if (optionalPlots_) {
2675  meRecoPosInSimCollection_->Fill(recopv.at(ir).sim);
2676  meRecoPosInRecoOrigCollection_->Fill(recopv.at(ir).OriginalIndex);
2677  }
2678  }
2679  if (simpv.at(is).eventId.bunchCrossing() == 0 && simpv.at(is).eventId.event() == 0) {
2680  if (!recopv.at(ir).is_signal()) {
2681  edm::LogWarning("Primary4DVertexValidation")
2682  << "Reco vtx leading match inconsistent: BX/ID " << simpv.at(is).eventId.bunchCrossing() << " "
2683  << simpv.at(is).eventId.event();
2684  }
2685  meRecoPVPosSignal_->Fill(
2686  recopv.at(ir).OriginalIndex); // position in reco vtx correction associated to sim signal
2687  if (!signal_is_highest_pt) {
2688  meRecoPVPosSignalNotHighestPt_->Fill(
2689  recopv.at(ir).OriginalIndex); // position in reco vtx correction associated to sim signal
2690  }
2691  }
2692  LogTrace("Primary4DVertexValidation") << "*** Matching RECO: " << ir << "with SIM: " << is << " ***";
2693  LogTrace("Primary4DVertexValidation") << "Match Quality is " << recopv.at(ir).matchQuality;
2694  LogTrace("Primary4DVertexValidation") << "****";
2695  }
2696  } // ndof
2697  }
2698  }
2699 
2700  // dz histos
2701  for (unsigned int iv = 0; iv < recVtxs->size() - 1; iv++) {
2702  if (recVtxs->at(iv).ndof() > selNdof_) {
2703  double mindistance_realreal = 1e10;
2704 
2705  for (unsigned int jv = iv; jv < recVtxs->size(); jv++) {
2706  if ((!(jv == iv)) && select(recVtxs->at(jv))) {
2707  double dz = recVtxs->at(iv).z() - recVtxs->at(jv).z();
2708  double dtsigma = std::sqrt(recVtxs->at(iv).covariance(3, 3) + recVtxs->at(jv).covariance(3, 3));
2709  double dt = (std::abs(dz) <= deltaZcut_ && dtsigma > 0.)
2710  ? (recVtxs->at(iv).t() - recVtxs->at(jv).t()) / dtsigma
2711  : -9999.;
2712  if (recopv.at(iv).is_real() && recopv.at(jv).is_real()) {
2713  meDeltaZrealreal_->Fill(std::abs(dz));
2714  if (dt != -9999.) {
2715  meDeltaTrealreal_->Fill(std::abs(dt));
2716  }
2717  if (std::abs(dz) < std::abs(mindistance_realreal)) {
2718  mindistance_realreal = dz;
2719  }
2720  } else if (recopv.at(iv).is_fake() && recopv.at(jv).is_fake()) {
2721  meDeltaZfakefake_->Fill(std::abs(dz));
2722  if (dt != -9999.) {
2723  meDeltaTfakefake_->Fill(std::abs(dt));
2724  }
2725  }
2726  }
2727  }
2728 
2729  double mindistance_fakereal = 1e10;
2730  double mindistance_realfake = 1e10;
2731  for (unsigned int jv = 0; jv < recVtxs->size(); jv++) {
2732  if ((!(jv == iv)) && select(recVtxs->at(jv))) {
2733  double dz = recVtxs->at(iv).z() - recVtxs->at(jv).z();
2734  double dtsigma = std::sqrt(recVtxs->at(iv).covariance(3, 3) + recVtxs->at(jv).covariance(3, 3));
2735  double dt = (std::abs(dz) <= deltaZcut_ && dtsigma > 0.)
2736  ? (recVtxs->at(iv).t() - recVtxs->at(jv).t()) / dtsigma
2737  : -9999.;
2738  if (recopv.at(iv).is_fake() && recopv.at(jv).is_real()) {
2739  meDeltaZfakereal_->Fill(std::abs(dz));
2740  if (dt != -9999.) {
2741  meDeltaTfakereal_->Fill(std::abs(dt));
2742  }
2743  if (std::abs(dz) < std::abs(mindistance_fakereal)) {
2744  mindistance_fakereal = dz;
2745  }
2746  }
2747 
2748  if (recopv.at(iv).is_real() && recopv.at(jv).is_fake() && (std::abs(dz) < std::abs(mindistance_realfake))) {
2749  mindistance_realfake = dz;
2750  }
2751  }
2752  }
2753  } // ndof
2754  }
2755 
2756 } // end of analyze
2757 
2758 void Primary4DVertexValidation::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
2759  edm::ParameterSetDescription desc;
2760 
2761  desc.add<std::string>("folder", "MTD/Vertices");
2762  desc.add<edm::InputTag>("TPtoRecoTrackAssoc", edm::InputTag("trackingParticleRecoTrackAsssociation"));
2763  desc.add<edm::InputTag>("mtdTracks", edm::InputTag("trackExtenderWithMTD"));
2764  desc.add<edm::InputTag>("SimTag", edm::InputTag("mix", "MergedTrackTruth"));
2765  desc.add<edm::InputTag>("offlineBS", edm::InputTag("offlineBeamSpot"));
2766  desc.add<edm::InputTag>("offline4DPV", edm::InputTag("offlinePrimaryVertices4D"));
2767  desc.add<edm::InputTag>("trackAssocSrc", edm::InputTag("trackExtenderWithMTD:generalTrackassoc"))
2768  ->setComment("Association between General and MTD Extended tracks");
2769  desc.add<edm::InputTag>("pathLengthSrc", edm::InputTag("trackExtenderWithMTD:generalTrackPathLength"));
2770  desc.add<edm::InputTag>("momentumSrc", edm::InputTag("trackExtenderWithMTD:generalTrackp"));
2771  desc.add<edm::InputTag>("tmtd", edm::InputTag("trackExtenderWithMTD:generalTracktmtd"));
2772  desc.add<edm::InputTag>("timeSrc", edm::InputTag("trackExtenderWithMTD:generalTracktmtd"));
2773  desc.add<edm::InputTag>("sigmaSrc", edm::InputTag("trackExtenderWithMTD:generalTracksigmatmtd"));
2774  desc.add<edm::InputTag>("t0PID", edm::InputTag("tofPID:t0"));
2775  desc.add<edm::InputTag>("sigmat0PID", edm::InputTag("tofPID:sigmat0"));
2776  desc.add<edm::InputTag>("t0SafePID", edm::InputTag("tofPID:t0safe"));
2777  desc.add<edm::InputTag>("sigmat0SafePID", edm::InputTag("tofPID:sigmat0safe"));
2778  desc.add<edm::InputTag>("trackMVAQual", edm::InputTag("mtdTrackQualityMVA:mtdQualMVA"));
2779  desc.add<edm::InputTag>("tofPi", edm::InputTag("trackExtenderWithMTD:generalTrackTofPi"));
2780  desc.add<edm::InputTag>("tofK", edm::InputTag("trackExtenderWithMTD:generalTrackTofK"));
2781  desc.add<edm::InputTag>("tofP", edm::InputTag("trackExtenderWithMTD:generalTrackTofP"));
2782  desc.add<edm::InputTag>("probPi", edm::InputTag("tofPID:probPi"));
2783  desc.add<edm::InputTag>("probK", edm::InputTag("tofPID:probK"));
2784  desc.add<edm::InputTag>("probP", edm::InputTag("tofPID:probP"));
2785  desc.add<bool>("useOnlyChargedTracks", true);
2786  desc.addUntracked<bool>("optionalPlots", false);
2787  desc.add<bool>("use3dNoTime", false);
2788  desc.add<double>("trackweightTh", 0.5);
2789  desc.add<double>("mvaTh", 0.8);
2790  desc.add<double>("minProbHeavy", 0.75);
2791 
2792  //lineDensity parameters have been obtained by fitting the distribution of the z position of the vertices,
2793  //using a 200k single mu ptGun sample (gaussian fit)
2794  std::vector<double> lDP;
2795  lDP.push_back(1.87);
2796  lDP.push_back(0.);
2797  lDP.push_back(42.5);
2798  desc.add<std::vector<double>>("lineDensityPar", lDP);
2799  descriptions.add("vertices4DValid", desc);
2800 }
2801 
2802 void Primary4DVertexValidation::printMatchedRecoTrackInfo(const reco::Vertex& vtx,
2803  const reco::TrackBaseRef& trk,
2804  const TrackingParticleRef& tp,
2805  const unsigned int& categ) {
2806  std::string strTrk;
2807  switch (categ) {
2808  case 0:
2809  strTrk = "Reco_Track:";
2810  break;
2811  case 1:
2812  strTrk = "SecRecoTrk:";
2813  break;
2814  case 2:
2815  strTrk = "PU_RecoTrk:";
2816  break;
2817  }
2818  LogTrace("Primary4DVertexValidation") << strTrk << " w =" << std::setw(6) << std::setprecision(2)
2819  << vtx.trackWeight(trk) << " pt =" << std::setw(6) << std::setprecision(2)
2820  << trk->pt() << " eta =" << std::setw(6) << std::setprecision(2) << trk->eta()
2821  << " MatchedTP: Pt =" << std::setw(6) << std::setprecision(2) << tp->pt()
2822  << " eta =" << std::setw(6) << std::setprecision(2) << tp->eta()
2823  << " Parent vtx: z =" << std::setw(8) << std::setprecision(4)
2824  << tp->parentVertex()->position().z() << " t =" << std::setw(8)
2825  << std::setprecision(4) << tp->parentVertex()->position().t() * simUnit_
2826  << " BX =" << tp->parentVertex()->eventId().bunchCrossing()
2827  << " ev =" << tp->parentVertex()->eventId().event() << std::endl;
2828 }
2829 
2830 void Primary4DVertexValidation::printSimVtxRecoVtxInfo(
2831  const struct Primary4DVertexValidation::simPrimaryVertex& simpVtx,
2832  const struct Primary4DVertexValidation::recoPrimaryVertex& recopVtx) {
2833  LogTrace("Primary4DVertexValidation") << "Sim vtx (x,y,z,t) = (" << std::setprecision(4) << simpVtx.x << ","
2834  << std::setprecision(4) << simpVtx.y << "," << std::setprecision(4) << simpVtx.z
2835  << "," << std::setprecision(4) << simpVtx.t * simUnit_ << ")"
2836  << " Simvtx.rec = " << simpVtx.rec;
2837  LogTrace("Primary4DVertexValidation") << "Sim vtx: pt = " << std::setprecision(4) << simpVtx.pt
2838  << " ptsq = " << std::setprecision(6) << simpVtx.ptsq
2839  << " nGenTrk = " << simpVtx.nGenTrk
2840  << " nmatch recotrks = " << simpVtx.num_matched_reco_tracks;
2841  LogTrace("Primary4DVertexValidation") << "Reco vtx (x,y,z) = (" << std::setprecision(4) << recopVtx.x << ","
2842  << std::setprecision(4) << recopVtx.y << "," << std::setprecision(4)
2843  << recopVtx.z << ")"
2844  << " Recovtx.sim = " << recopVtx.sim;
2845  LogTrace("Primary4DVertexValidation") << "Reco vtx: pt = " << std::setprecision(4) << recopVtx.pt
2846  << " ptsq = " << std::setprecision(6) << recopVtx.ptsq
2847  << " nrecotrks = " << recopVtx.nRecoTrk
2848  << " nmatch simtrks = " << recopVtx.num_matched_sim_tracks;
2849  LogTrace("Primary4DVertexValidation") << "wnt " << recopVtx.sumwnt << " wos = " << recopVtx.sumwos;
2850  for (auto iTP = simpVtx.sim_vertex->daughterTracks_begin(); iTP != simpVtx.sim_vertex->daughterTracks_end(); ++iTP) {
2851  if (use_only_charged_tracks_ && (**iTP).charge() == 0) {
2852  continue;
2853  }
2854  LogTrace("Primary4DVertexValidation")
2855  << "Daughter track of sim vertex: pt =" << std::setw(6) << std::setprecision(2) << (*iTP)->pt()
2856  << " eta =" << std::setw(6) << std::setprecision(2) << (*iTP)->eta();
2857  }
2858 }
2859 
2860 const bool Primary4DVertexValidation::trkTPSelLV(const TrackingParticle& tp) {
2861  bool match = false;
2862  if (tp.status() != 1) {
2863  return match;
2864  }
2865 
2866  auto x_pv = tp.parentVertex()->position().x();
2867  auto y_pv = tp.parentVertex()->position().y();
2868  auto z_pv = tp.parentVertex()->position().z();
2869 
2870  auto r_pv = std::sqrt(x_pv * x_pv + y_pv * y_pv);
2871 
2872  match =
2873  tp.charge() != 0 && std::abs(tp.eta()) < etacutGEN_ && tp.pt() > pTcut_ && r_pv < rBTL_ && std::abs(z_pv) < zETL_;
2874  return match;
2875 }
2876 
2877 const bool Primary4DVertexValidation::trkRecSel(const reco::TrackBase& trk) {
2878  bool match = false;
2879  match = std::abs(trk.eta()) <= etacutREC_ && trk.pt() > pTcut_;
2880  return match;
2881 }
2882 
2883 DEFINE_FWK_MODULE(Primary4DVertexValidation);
Primary4DVertexValidation::trackMaxBtlEta_
static constexpr double trackMaxBtlEta_
Definition: Primary4DVertexValidation.cc:274
dt
float dt
Definition: AMPTWrapper.h:136
Primary4DVertexValidation::Rec4DVerToken_
edm::EDGetTokenT< edm::View< reco::Vertex > > Rec4DVerToken_
Definition: Primary4DVertexValidation.cc:309
ProducerED_cfi.InputTag
InputTag
Definition: ProducerED_cfi.py:5
deDxTools::esConsumes
ESGetTokenH3DDVariant esConsumes(std::string const &Record, edm::ConsumesCollector &)
Definition: DeDxTools.cc:283
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
Definition: ParameterSet.h:307
Primary4DVertexValidation::trackingVertexCollectionToken_
edm::EDGetTokenT< TrackingVertexCollection > trackingVertexCollectionToken_
Definition: Primary4DVertexValidation.cc:305
bphysicsOniaDQM_cfi.vertex
vertex
Definition: bphysicsOniaDQM_cfi.py:7
Primary4DVertexValidation::trackMaxEtlEta_
static constexpr double trackMaxEtlEta_
Definition: Primary4DVertexValidation.cc:276
genfragment_ptgun_cfg.PartID
PartID
Definition: genfragment_ptgun_cfg.py:6
mps_fire.i
i
Definition: mps_fire.py:429
Primary4DVertexValidation::etacutREC_
static constexpr double etacutREC_
Definition: Primary4DVertexValidation.cc:271
l1tEGammaCrystalsEmulatorProducer_cfi.scale
scale
Definition: l1tEGammaCrystalsEmulatorProducer_cfi.py:10
folder_
std::string folder_
Definition: DQMEDAnalyzer.cc:60
muonTagProbeFilters_cff.matched
matched
Definition: muonTagProbeFilters_cff.py:62
Primary4DVertexValidation::mePUTrackRecLVRelSumWntvsSumWnt_
MonitorElement * mePUTrackRecLVRelSumWntvsSumWnt_
Definition: Primary4DVertexValidation.cc:426
Primary4DVertexValidation::trackingParticleCollectionToken_
edm::EDGetTokenT< TrackingParticleCollection > trackingParticleCollectionToken_
Definition: Primary4DVertexValidation.cc:304
dqm::implementation::NavigatorBase::setCurrentFolder
virtual void setCurrentFolder(std::string const &fullpath)
Definition: DQMStore.cc:36
Primary4DVertexValidation::sigmat0SafePidToken_
edm::EDGetTokenT< edm::ValueMap< float > > sigmat0SafePidToken_
Definition: Primary4DVertexValidation.cc:319
ParticleDataTable
HepPDT::ParticleDataTable ParticleDataTable
Definition: ParticleDataTable.h:8
Primary4DVertexValidation::getSimPVs
std::vector< Primary4DVertexValidation::simPrimaryVertex > getSimPVs(const edm::Handle< TrackingVertexCollection > &)
Definition: Primary4DVertexValidation.cc:1586
Primary4DVertexValidation::RecBeamSpotToken_
edm::EDGetTokenT< reco::BeamSpot > RecBeamSpotToken_
Definition: Primary4DVertexValidation.cc:308
edm::isNotFinite
constexpr bool isNotFinite(T x)
Definition: isFinite.h:9
Primary4DVertexValidation::trkTPSelLV
const bool trkTPSelLV(const TrackingParticle &)
Definition: Primary4DVertexValidation.cc:2860
Primary4DVertexValidation::trkRecSel
const bool trkRecSel(const reco::TrackBase &)
Definition: Primary4DVertexValidation.cc:2877
Primary4DVertexValidation::select
bool select(const reco::Vertex &, int level=0)
Definition: Primary4DVertexValidation.cc:1438
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:1752
Primary4DVertexValidation::observablesFromJets
void observablesFromJets(const std::vector< reco::Track > &, const std::vector< double > &, const std::vector< int > &, const std::string &, unsigned int &, double &, double &, double &, double &)
Definition: Primary4DVertexValidation.cc:1457
Primary4DVertexValidation::probKToken_
edm::EDGetTokenT< edm::ValueMap< float > > probKToken_
Definition: Primary4DVertexValidation.cc:327
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
reco::VertexBaseRef
edm::RefToBase< reco::Vertex > VertexBaseRef
persistent reference to a Vertex, using views
Definition: VertexFwd.h:21
edm::Handle::product
T const * product() const
Definition: Handle.h:70
L1BJetProducer_cff.vtx
vtx
Definition: L1BJetProducer_cff.py:10
Primary4DVertexValidation::matchRecoTrack2SimSignal
bool matchRecoTrack2SimSignal(const reco::TrackBaseRef &)
Definition: Primary4DVertexValidation.cc:1380
caHitNtupletGeneratorKernels::iev
TkSoAView< TrackerTraits > HitToTuple< TrackerTraits > const *__restrict__ int32_t int32_t int iev
Definition: CAHitNtupletGeneratorKernelsImpl.h:860
Primary4DVertexValidation::etacutGEN_
static constexpr double etacutGEN_
Definition: Primary4DVertexValidation.cc:270
Primary4DVertexValidation::sigmat0PidToken_
edm::EDGetTokenT< edm::ValueMap< float > > sigmat0PidToken_
Definition: Primary4DVertexValidation.cc:317
l1ctLayer1_cff.dZ
dZ
Definition: l1ctLayer1_cff.py:69
PDWG_EXODelayedJetMET_cff.jets
jets
Definition: PDWG_EXODelayedJetMET_cff.py:14
Primary4DVertexValidation::minThrSumWos_
static constexpr double minThrSumWos_
Definition: Primary4DVertexValidation.cc:279
Primary4DVertexValidation::meFakeTrackRecLVRelSumWosvsSumWos_
MonitorElement * meFakeTrackRecLVRelSumWosvsSumWos_
Definition: Primary4DVertexValidation.cc:429
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
pfDeepBoostedJetPreprocessParams_cfi.sv
sv
Definition: pfDeepBoostedJetPreprocessParams_cfi.py:352
Primary4DVertexValidation::meJetsRecLVPURelSumPt2vsSumPt2_
MonitorElement * meJetsRecLVPURelSumPt2vsSumPt2_
Definition: Primary4DVertexValidation.cc:466
ALPAKA_ACCELERATOR_NAMESPACE::brokenline::constexpr
if constexpr(n > 3)
Definition: BrokenLine.h:164
nJets
static constexpr int nJets
Definition: Phase2L1CaloJetEmulator.h:16
cms::cuda::assert
assert(be >=bs)
Primary4DVertexValidation::analyze
void analyze(const edm::Event &, const edm::EventSetup &) override
Definition: Primary4DVertexValidation.cc:1976
Primary4DVertexValidation::timeToken_
edm::EDGetTokenT< edm::ValueMap< float > > timeToken_
Definition: Primary4DVertexValidation.cc:314
Primary4DVertexValidation::getRecoPVs
std::vector< Primary4DVertexValidation::recoPrimaryVertex > getRecoPVs(const edm::Handle< edm::View< reco::Vertex >> &)
Definition: Primary4DVertexValidation.cc:1688
hcalRecHitTable_cff.time
time
Definition: hcalRecHitTable_cff.py:14
Primary4DVertexValidation::trackMinEtlEta_
static constexpr double trackMinEtlEta_
Definition: Primary4DVertexValidation.cc:275
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
edm::AssociationMap::find
const_iterator find(const key_type &k) const
find element with specified reference key
Definition: AssociationMap.h:175
LHEGenericFilter_cfi.ParticleID
ParticleID
Definition: LHEGenericFilter_cfi.py:6
LogTrace
#define LogTrace(id)
Definition: MessageLogger.h:242
Primary4DVertexValidation::meTrackResMassProtons_
MonitorElement * meTrackResMassProtons_[3]
Definition: Primary4DVertexValidation.cc:482
Primary4DVertexValidation::meTrackResMassTrueProtons_
MonitorElement * meTrackResMassTrueProtons_[3]
Definition: Primary4DVertexValidation.cc:483
edm::AssociationMap::end
const_iterator end() const
last iterator over the map (read only)
Definition: AssociationMap.h:173
dqm::impl::MonitorElement::Fill
void Fill(long long x)
Definition: MonitorElement.h:292
edm::second
U second(std::pair< T, U > const &p)
Definition: ParameterSet.cc:222
math::XYZTLorentzVector
XYZTLorentzVectorD XYZTLorentzVector
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:29
reco::TrackBase::pt
double pt() const
track transverse momentum
Definition: TrackBase.h:637
Primary4DVertexValidation::timeFromTrueMass
double timeFromTrueMass(double, double, double, double)
Definition: Primary4DVertexValidation.cc:1426
iEvent
int iEvent
Definition: GenABIO.cc:224
Primary4DVertexValidation::LorentzVector
math::XYZTLorentzVector LorentzVector
Definition: Primary4DVertexValidation.cc:67
ALPAKA_ACCELERATOR_NAMESPACE::vertexFinder::iv
int32_t *__restrict__ iv
Definition: clusterTracksByDensity.h:58
Primary4DVertexValidation::meSecTrackRecLVRelSumWosvsSumWos_
MonitorElement * meSecTrackRecLVRelSumWosvsSumWos_
Definition: Primary4DVertexValidation.cc:428
Primary4DVertexValidation::simPrimaryVertex::addTrack
void addTrack(unsigned int irecv, double twos, double twt)
Definition: Primary4DVertexValidation.cc:113
dqm::implementation::IBooker::bookProfile
MonitorElement * bookProfile(TString const &name, TString const &title, int nchX, double lowX, double highX, int, double lowY, double highY, char const *option="s", FUNC onbooking=NOOP())
Definition: DQMStore.h:408
reco::TrackBase::dzError
double dzError() const
error on dz
Definition: TrackBase.h:778
V0Monitor_cff.v0
v0
Definition: V0Monitor_cff.py:7
Primary4DVertexValidation::fillDescriptions
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
Definition: Primary4DVertexValidation.cc:2758
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:23
Primary4DVertexValidation::Primary4DVertexValidation
Primary4DVertexValidation(const edm::ParameterSet &)
Definition: Primary4DVertexValidation.cc:523
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
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:97
Primary4DVertexValidation::probPiToken_
edm::EDGetTokenT< edm::ValueMap< float > > probPiToken_
Definition: Primary4DVertexValidation.cc:326
Primary4DVertexValidation::minThrMetPt_
static constexpr double minThrMetPt_
Definition: Primary4DVertexValidation.cc:282
funct::tan
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
Primary4DVertexValidation::mePUTrackRecLVRelSumPt2vsSumPt2_
MonitorElement * mePUTrackRecLVRelSumPt2vsSumPt2_
Definition: Primary4DVertexValidation.cc:431
Primary4DVertexValidation::probPToken_
edm::EDGetTokenT< edm::ValueMap< float > > probPToken_
Definition: Primary4DVertexValidation.cc:328
Primary4DVertexValidation::meJetsRecLVFakeRelSumPt2vsSumPt2_
MonitorElement * meJetsRecLVFakeRelSumPt2vsSumPt2_
Definition: Primary4DVertexValidation.cc:467
reco::TrackBase::phi
double phi() const
azimuthal angle of momentum vector
Definition: TrackBase.h:649
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
Primary4DVertexValidation::momentumToken_
edm::EDGetTokenT< edm::ValueMap< float > > momentumToken_
Definition: Primary4DVertexValidation.cc:313
Primary4DVertexValidation::meTrackResMassTruePions_
MonitorElement * meTrackResMassTruePions_[3]
Definition: Primary4DVertexValidation.cc:485
Primary4DVertexValidation::tofKToken_
edm::EDGetTokenT< edm::ValueMap< float > > tofKToken_
Definition: Primary4DVertexValidation.cc:324
Primary4DVertexValidation::RecTrackToken_
edm::EDGetTokenT< reco::TrackCollection > RecTrackToken_
Definition: Primary4DVertexValidation.cc:300
HLT_2024v14_cff.distance
distance
Definition: HLT_2024v14_cff.py:5740
Primary4DVertexValidation::lineDensityPar_
const std::vector< double > lineDensityPar_
Definition: Primary4DVertexValidation.cc:296
ParticleData
HepPDT::ParticleData ParticleData
Definition: ParticleDataTable.h:9
edm::EventSetup::getHandle
ESHandle< T > getHandle(const ESGetToken< T, R > &iToken) const
Definition: EventSetup.h:130
Primary4DVertexValidation::r2s_
const reco::RecoToSimCollection * r2s_
Definition: Primary4DVertexValidation.cc:297
Primary4DVertexValidation::vecPileupSummaryInfoToken_
edm::EDGetTokenT< std::vector< PileupSummaryInfo > > vecPileupSummaryInfoToken_
Definition: Primary4DVertexValidation.cc:302
nt
int nt
Definition: AMPTWrapper.h:42
Primary4DVertexValidation::simPrimaryVertex::simPrimaryVertex
simPrimaryVertex(double x1, double y1, double z1, double t1)
Definition: Primary4DVertexValidation.cc:71
Primary4DVertexValidation::deltaZcut_
static constexpr double deltaZcut_
Definition: Primary4DVertexValidation.cc:273
EgHLTOffHistBins_cfi.mass
mass
Definition: EgHLTOffHistBins_cfi.py:34
Primary4DVertexValidation::t0PidToken_
edm::EDGetTokenT< edm::ValueMap< float > > t0PidToken_
Definition: Primary4DVertexValidation.cc:316
TrackingVertexRef
edm::Ref< TrackingVertexCollection > TrackingVertexRef
Definition: TrackingVertexContainer.h:9
Primary4DVertexValidation::NOT_MATCHED
static constexpr unsigned int NOT_MATCHED
Definition: Primary4DVertexValidation.cc:261
reco::TrackBase::eta
double eta() const
pseudorapidity of momentum vector
Definition: TrackBase.h:652
Primary4DVertexValidation::simToRecoAssociationToken_
edm::EDGetTokenT< reco::SimToRecoCollection > simToRecoAssociationToken_
Definition: Primary4DVertexValidation.cc:306
Primary4DVertexValidation::trackAssocToken_
edm::EDGetTokenT< edm::ValueMap< int > > trackAssocToken_
Definition: Primary4DVertexValidation.cc:311
angle_units::operators::convertMmToCm
constexpr NumType convertMmToCm(NumType millimeters)
Definition: angle_units.h:44
Primary4DVertexValidation::minThrSumPz_
static constexpr double minThrSumPz_
Definition: Primary4DVertexValidation.cc:283
Primary4DVertexValidation::isParticle
void isParticle(const reco::TrackBaseRef &, const edm::ValueMap< float > &, const edm::ValueMap< float > &, const edm::ValueMap< float > &, const edm::ValueMap< float > &, const edm::ValueMap< float > &, unsigned int &, bool &, bool &, bool &, bool &)
Definition: Primary4DVertexValidation.cc:1530
edm::ConfigurationDescriptions::add
void add(std::string const &label, ParameterSetDescription const &psetDescription)
Definition: ConfigurationDescriptions.cc:57
nano_mu_local_reco_cff.bool
bool
Definition: nano_mu_local_reco_cff.py:14
Primary4DVertexValidation::minThrSumPt2_
static constexpr double minThrSumPt2_
Definition: Primary4DVertexValidation.cc:281
testProducerWithPsetDescEmpty_cfi.y1
y1
Definition: testProducerWithPsetDescEmpty_cfi.py:29
Primary4DVertexValidation::recoPrimaryVertex::addTrack
void addTrack(unsigned int iev, double twos, double wt)
Definition: Primary4DVertexValidation.cc:183
edm::HandleBase::isValid
bool isValid() const
Definition: HandleBase.h:70
Primary4DVertexValidation::printSimVtxRecoVtxInfo
void printSimVtxRecoVtxInfo(const struct Primary4DVertexValidation::simPrimaryVertex &, const struct Primary4DVertexValidation::recoPrimaryVertex &)
Definition: Primary4DVertexValidation.cc:2830
BTaggingMonitor_cfi.met
met
Definition: BTaggingMonitor_cfi.py:74
reco::TrackBase::theta
double theta() const
polar angle
Definition: TrackBase.h:602
Primary4DVertexValidation::getMatchedTP
std::pair< const edm::Ref< std::vector< TrackingParticle > > *, int > getMatchedTP(const reco::TrackBaseRef &, const TrackingVertexRef &)
Definition: Primary4DVertexValidation.cc:1397
reco
fixed size matrix
Definition: AlignmentAlgorithmBase.h:46
edm
HLT enums.
Definition: AlignableModifier.h:19
Primary4DVertexValidation::s2r_
const reco::SimToRecoCollection * s2r_
Definition: Primary4DVertexValidation.cc:298
Primary4DVertexValidation::tofPiToken_
edm::EDGetTokenT< edm::ValueMap< float > > tofPiToken_
Definition: Primary4DVertexValidation.cc:323
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:567
Primary4DVertexValidation::recoToSimAssociationToken_
edm::EDGetTokenT< reco::RecoToSimCollection > recoToSimAssociationToken_
Definition: Primary4DVertexValidation.cc:307
match
std::pair< typename Association::data_type::first_type, double > match(Reference key, Association association, bool bestMatchByMaxValue)
Generic matching function.
Definition: Utils.h:10
gather_cfg.cout
cout
Definition: gather_cfg.py:144
Primary4DVertexValidation::pdtToken_
edm::ESGetToken< HepPDT::ParticleDataTable, edm::DefaultRecord > pdtToken_
Definition: Primary4DVertexValidation.cc:329
JetChargeProducer_cfi.exp
exp
Definition: JetChargeProducer_cfi.py:6
edm::LogWarning
Log< level::Warning, false > LogWarning
Definition: MessageLogger.h:128
dqm::implementation::IBooker::book1D
MonitorElement * book1D(TString const &name, TString const &title, int const nchX, double const lowX, double const highX, FUNC onbooking=NOOP())
Definition: DQMStore.h:98
Primary4DVertexValidation::t0SafePidToken_
edm::EDGetTokenT< edm::ValueMap< float > > t0SafePidToken_
Definition: Primary4DVertexValidation.cc:318
Primary4DVertexValidation::tmtdToken_
edm::EDGetTokenT< edm::ValueMap< float > > tmtdToken_
Definition: Primary4DVertexValidation.cc:322
testProducerWithPsetDescEmpty_cfi.x1
x1
Definition: testProducerWithPsetDescEmpty_cfi.py:33
Primary4DVertexValidation::zWosMatchMax_
static constexpr double zWosMatchMax_
Definition: Primary4DVertexValidation.cc:269
Primary4DVertexValidation::mePUTrackRecLVRelSumWosvsSumWos_
MonitorElement * mePUTrackRecLVRelSumWosvsSumWos_
Definition: Primary4DVertexValidation.cc:427
Primary4DVertexValidation::minThrSumPt_
static constexpr double minThrSumPt_
Definition: Primary4DVertexValidation.cc:280
Primary4DVertexValidation::tofPToken_
edm::EDGetTokenT< edm::ValueMap< float > > tofPToken_
Definition: Primary4DVertexValidation.cc:325
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
TtFullHadEvtBuilder_cfi.sumPt
sumPt
Definition: TtFullHadEvtBuilder_cfi.py:38
Primary4DVertexValidation::minThrSumWnt_
static constexpr double minThrSumWnt_
Definition: Primary4DVertexValidation.cc:278
Primary4DVertexValidation::getWosWnt
void getWosWnt(const reco::Vertex &, const reco::TrackBaseRef &, const edm::ValueMap< float > &, const edm::ValueMap< float > &, const edm::Handle< reco::BeamSpot > &, double &, double &)
Definition: Primary4DVertexValidation.cc:1560
edm::Run
Definition: Run.h:45
Primary4DVertexValidation::trackMVAQualToken_
edm::EDGetTokenT< edm::ValueMap< float > > trackMVAQualToken_
Definition: Primary4DVertexValidation.cc:320
Primary4DVertexValidation::printMatchedRecoTrackInfo
void printMatchedRecoTrackInfo(const reco::Vertex &, const reco::TrackBaseRef &, const TrackingParticleRef &, const unsigned int &)
Definition: Primary4DVertexValidation.cc:2802
Primary4DVertexValidation::pathLengthToken_
edm::EDGetTokenT< edm::ValueMap< float > > pathLengthToken_
Definition: Primary4DVertexValidation.cc:312
Generated for CMSSW Reference Manual by doxygen 1.8.14