Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-01-29 23:45:29 

0001 #include <numeric>
0002 #include <string>
0003 #include <vector>
0004 #include <map>
0005 #include <algorithm>
0006 
0007 #include "FWCore/Framework/interface/Frameworkfwd.h"
0008 #include "FWCore/Framework/interface/ESHandle.h"
0009 
0010 #include "DataFormats/Common/interface/ValidHandle.h"
0011 #include "FWCore/Framework/interface/Event.h"
0012 #include "FWCore/Framework/interface/EventSetup.h"
0013 #include "FWCore/Framework/interface/MakerMacros.h"
0014 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0015 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0016 
0017 #include "DataFormats/Math/interface/LorentzVector.h"
0018 #include "DataFormats/Math/interface/Point3D.h"
0019 
0020 // reco track and vertex
0021 #include "DataFormats/TrackReco/interface/Track.h"
0022 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0023 #include "DataFormats/VertexReco/interface/Vertex.h"
0024 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0025 #include "RecoVertex/VertexPrimitives/interface/TransientVertex.h"
0026 
0027 // TrackingParticle
0028 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
0029 #include "SimDataFormats/Associations/interface/TrackToTrackingParticleAssociator.h"
0030 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h"
0031 #include "SimDataFormats/TrackingAnalysis/interface/TrackingVertexContainer.h"
0032 
0033 // pile-up
0034 #include "SimDataFormats/PileupSummaryInfo/interface/PileupSummaryInfo.h"
0035 
0036 // associator
0037 #include "SimTracker/VertexAssociation/interface/calculateVertexSharedTracks.h"
0038 
0039 // vertexing
0040 #include "RecoVertex/PrimaryVertexProducer/interface/TrackFilterForPVFinding.h"
0041 
0042 // simulated vertex
0043 #include "SimDataFormats/Associations/interface/VertexToTrackingVertexAssociator.h"
0044 
0045 // DQM
0046 #include "DQMServices/Core/interface/DQMEDAnalyzer.h"
0047 #include "DQMServices/Core/interface/DQMStore.h"
0048 
0049 #include "DataFormats/Math/interface/deltaR.h"
0050 #include "DataFormats/Math/interface/GeantUnits.h"
0051 #include "CLHEP/Units/PhysicalConstants.h"
0052 
0053 //class declaration
0054 class Primary4DVertexValidation : public DQMEDAnalyzer {
0055   typedef math::XYZTLorentzVector LorentzVector;
0056 
0057   // auxiliary class holding simulated vertices
0058   struct simPrimaryVertex {
0059     simPrimaryVertex(double x1, double y1, double z1, double t1)
0060         : x(x1),
0061           y(y1),
0062           z(z1),
0063           t(t1),
0064           ptsq(0),
0065           closest_vertex_distance_z(-1.),
0066           nGenTrk(0),
0067           num_matched_reco_tracks(0),
0068           average_match_quality(0.0) {
0069       ptot.setPx(0);
0070       ptot.setPy(0);
0071       ptot.setPz(0);
0072       ptot.setE(0);
0073       p4 = LorentzVector(0, 0, 0, 0);
0074       r = sqrt(x * x + y * y);
0075     };
0076     double x, y, z, r, t;
0077     HepMC::FourVector ptot;
0078     LorentzVector p4;
0079     double ptsq;
0080     double closest_vertex_distance_z;
0081     int nGenTrk;
0082     int num_matched_reco_tracks;
0083     float average_match_quality;
0084     EncodedEventId eventId;
0085     TrackingVertexRef sim_vertex;
0086     int OriginalIndex = -1;
0087 
0088     unsigned int nwosmatch = 0;                    // number of recvertices dominated by this simevt (by wos)
0089     unsigned int nwntmatch = 0;                    // number of recvertices dominated by this simevt  (by nt)
0090     std::vector<unsigned int> wos_dominated_recv;  // list of dominated recv (by wos, size==nwosmatch)
0091 
0092     std::map<unsigned int, double> wnt;  // weighted number of tracks in recvtx (by index)
0093     std::map<unsigned int, double> wos;  // sum of wos in recvtx (by index) // oops -> this was int before 04-22
0094     double sumwos = 0;                   // sum of wos in any recvtx
0095     double sumwnt = 0;                   // sum of weighted tracks
0096     unsigned int rec = NOT_MATCHED;      // best match  (NO_MATCH if not matched)
0097     unsigned int matchQuality = 0;       // quality flag
0098 
0099     void addTrack(unsigned int irecv, double twos, double twt) {
0100       sumwnt += twt;
0101       if (wnt.find(irecv) == wnt.end()) {
0102         wnt[irecv] = twt;
0103       } else {
0104         wnt[irecv] += twt;
0105       }
0106 
0107       sumwos += twos;
0108       if (wos.find(irecv) == wos.end()) {
0109         wos[irecv] = twos;
0110       } else {
0111         wos[irecv] += twos;
0112       }
0113     }
0114   };
0115 
0116   // auxiliary class holding reconstructed vertices
0117   struct recoPrimaryVertex {
0118     recoPrimaryVertex(double x1, double y1, double z1)
0119         : x(x1),
0120           y(y1),
0121           z(z1),
0122           pt(0),
0123           ptsq(0),
0124           closest_vertex_distance_z(-1.),
0125           nRecoTrk(0),
0126           num_matched_sim_tracks(0),
0127           recVtx(nullptr) {
0128       r = sqrt(x * x + y * y);
0129     };
0130     double x, y, z, r;
0131     double pt;
0132     double ptsq;
0133     double closest_vertex_distance_z;
0134     int nRecoTrk;
0135     int num_matched_sim_tracks;
0136     const reco::Vertex* recVtx;
0137     reco::VertexBaseRef recVtxRef;
0138     int OriginalIndex = -1;
0139 
0140     std::map<unsigned int, double> wos;  // simevent -> wos
0141     std::map<unsigned int, double> wnt;  // simevent -> weighted number of truth matched tracks
0142     unsigned int wosmatch;               // index of the simevent providing the largest contribution to wos
0143     unsigned int wntmatch;               // index of the simevent providing the highest number of tracks
0144     double sumwos = 0;                   // total sum of wos of all truth matched tracks
0145     double sumwnt = 0;                   // total weighted number of truth matchted tracks
0146     double maxwos = 0;                   // largest wos sum from one sim event (wosmatch)
0147     double maxwnt = 0;                   // largest weighted  number of tracks from one sim event (ntmatch)
0148     int maxwosnt = 0;                    // number of tracks from the simevt with highest wos
0149     unsigned int sim = NOT_MATCHED;      // best match  (NO_MATCH if not matched)
0150     unsigned int matchQuality = 0;       // quality flag
0151 
0152     bool is_real() { return (matchQuality > 0) && (matchQuality < 99); }
0153 
0154     bool is_fake() { return (matchQuality <= 0) || (matchQuality >= 99); }
0155 
0156     bool is_signal() { return (sim == 0); }
0157 
0158     int split_from() {
0159       if (is_real())
0160         return -1;
0161       if ((maxwos > 0) && (maxwos > 0.3 * sumwos))
0162         return wosmatch;
0163       return -1;
0164     }
0165     bool other_fake() { return (is_fake() && (split_from() < 0)); }
0166 
0167     void addTrack(unsigned int iev, double twos, double wt) {
0168       sumwnt += wt;
0169       if (wnt.find(iev) == wnt.end()) {
0170         wnt[iev] = wt;
0171       } else {
0172         wnt[iev] += wt;
0173       }
0174 
0175       sumwos += twos;
0176       if (wos.find(iev) == wos.end()) {
0177         wos[iev] = twos;
0178       } else {
0179         wos[iev] += twos;
0180       }
0181     }
0182   };
0183 
0184 public:
0185   explicit Primary4DVertexValidation(const edm::ParameterSet&);
0186   ~Primary4DVertexValidation() override;
0187 
0188   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0189   void analyze(const edm::Event&, const edm::EventSetup&) override;
0190   void bookHistograms(DQMStore::IBooker& i, edm::Run const&, edm::EventSetup const&) override;
0191 
0192 private:
0193   void matchReco2Sim(std::vector<recoPrimaryVertex>&,
0194                      std::vector<simPrimaryVertex>&,
0195                      const edm::ValueMap<float>&,
0196                      const edm::ValueMap<float>&,
0197                      const edm::Handle<reco::BeamSpot>&);
0198   bool matchRecoTrack2SimSignal(const reco::TrackBaseRef&);
0199   const edm::Ref<std::vector<TrackingParticle>>* getMatchedTP(const reco::TrackBaseRef&, const TrackingVertexRef&);
0200   double timeFromTrueMass(double, double, double, double);
0201   bool select(const reco::Vertex&, int level = 0);
0202   std::vector<Primary4DVertexValidation::simPrimaryVertex> getSimPVs(const edm::Handle<TrackingVertexCollection>&);
0203   std::vector<Primary4DVertexValidation::recoPrimaryVertex> getRecoPVs(const edm::Handle<edm::View<reco::Vertex>>&);
0204 
0205   const bool mvaTPSel(const TrackingParticle&);
0206   const bool mvaRecSel(const reco::TrackBase&, const reco::Vertex&, const double&, const double&);
0207 
0208   // ----------member data ---------------------------
0209 
0210   const std::string folder_;
0211   static constexpr unsigned int NOT_MATCHED = 66666;
0212   static constexpr double simUnit_ = 1e9;     //sim time in s while reco time in ns
0213   static constexpr double c_ = 2.99792458e1;  //c in cm/ns
0214   static constexpr double mvaL_ = 0.5;        //MVA cuts for MVA categories
0215   static constexpr double mvaH_ = 0.8;
0216   static constexpr double selNdof_ = 4.;
0217   static constexpr double maxRank_ = 8.;
0218   static constexpr double maxTry_ = 10.;
0219   static constexpr double zWosMatchMax_ = 1.;
0220   static constexpr double etacutGEN_ = 4.;   // |eta| < 4;
0221   static constexpr double etacutREC_ = 3.;   // |eta| < 3;
0222   static constexpr double pTcut_ = 0.7;      // PT > 0.7 GeV
0223   static constexpr double deltaZcut_ = 0.1;  // dz separation 1 mm
0224   static constexpr double trackMaxBtlEta_ = 1.5;
0225   static constexpr double trackMinEtlEta_ = 1.6;
0226   static constexpr double trackMaxEtlEta_ = 3.;
0227   static constexpr double tol_ = 1.e-4;  // tolerance on reconstructed track time, [ns]
0228 
0229   static constexpr float c_cm_ns = geant_units::operators::convertMmToCm(CLHEP::c_light);  // [mm/ns] -> [cm/ns]
0230 
0231   bool use_only_charged_tracks_;
0232   bool debug_;
0233   bool optionalPlots_;
0234 
0235   const double minProbHeavy_;
0236   const double trackweightTh_;
0237   const double mvaTh_;
0238   const std::vector<double> lineDensityPar_;
0239   const reco::RecoToSimCollection* r2s_;
0240   const reco::SimToRecoCollection* s2r_;
0241 
0242   edm::EDGetTokenT<reco::TrackCollection> RecTrackToken_;
0243 
0244   edm::EDGetTokenT<std::vector<PileupSummaryInfo>> vecPileupSummaryInfoToken_;
0245 
0246   edm::EDGetTokenT<TrackingParticleCollection> trackingParticleCollectionToken_;
0247   edm::EDGetTokenT<TrackingVertexCollection> trackingVertexCollectionToken_;
0248   edm::EDGetTokenT<reco::SimToRecoCollection> simToRecoAssociationToken_;
0249   edm::EDGetTokenT<reco::RecoToSimCollection> recoToSimAssociationToken_;
0250   edm::EDGetTokenT<reco::BeamSpot> RecBeamSpotToken_;
0251   edm::EDGetTokenT<edm::View<reco::Vertex>> Rec4DVerToken_;
0252 
0253   edm::EDGetTokenT<edm::ValueMap<int>> trackAssocToken_;
0254   edm::EDGetTokenT<edm::ValueMap<float>> pathLengthToken_;
0255   edm::EDGetTokenT<edm::ValueMap<float>> momentumToken_;
0256   edm::EDGetTokenT<edm::ValueMap<float>> timeToken_;
0257 
0258   edm::EDGetTokenT<edm::ValueMap<float>> t0PidToken_;
0259   edm::EDGetTokenT<edm::ValueMap<float>> t0SafePidToken_;
0260   edm::EDGetTokenT<edm::ValueMap<float>> sigmat0SafePidToken_;
0261   edm::EDGetTokenT<edm::ValueMap<float>> trackMVAQualToken_;
0262 
0263   edm::EDGetTokenT<edm::ValueMap<float>> tmtdToken_;
0264   edm::EDGetTokenT<edm::ValueMap<float>> tofPiToken_;
0265   edm::EDGetTokenT<edm::ValueMap<float>> tofKToken_;
0266   edm::EDGetTokenT<edm::ValueMap<float>> tofPToken_;
0267   edm::EDGetTokenT<edm::ValueMap<float>> probPiToken_;
0268   edm::EDGetTokenT<edm::ValueMap<float>> probKToken_;
0269   edm::EDGetTokenT<edm::ValueMap<float>> probPToken_;
0270 
0271   //histogram declaration
0272   MonitorElement* meMVATrackEffPtTot_;
0273   MonitorElement* meMVATrackMatchedEffPtTot_;
0274   MonitorElement* meMVATrackMatchedEffPtMtd_;
0275   MonitorElement* meMVATrackEffEtaTot_;
0276   MonitorElement* meMVATrackMatchedEffEtaTot_;
0277   MonitorElement* meMVATrackMatchedEffEtaMtd_;
0278   MonitorElement* meMVATrackResTot_;
0279   MonitorElement* meMVATrackPullTot_;
0280   MonitorElement* meTrackResTot_;
0281   MonitorElement* meTrackPullTot_;
0282   MonitorElement* meTrackRes_[3];
0283   MonitorElement* meTrackPull_[3];
0284   MonitorElement* meTrackResMass_[3];
0285   MonitorElement* meTrackResMassTrue_[3];
0286   MonitorElement* meMVATrackZposResTot_;
0287   MonitorElement* meTrackZposResTot_;
0288   MonitorElement* meTrackZposRes_[3];
0289   MonitorElement* meTrack3DposRes_[3];
0290   MonitorElement* meTimeRes_;
0291   MonitorElement* meTimePull_;
0292   MonitorElement* meTimeSignalRes_;
0293   MonitorElement* meTimeSignalPull_;
0294   MonitorElement* mePUvsRealV_;
0295   MonitorElement* mePUvsOtherFakeV_;
0296   MonitorElement* mePUvsSplitV_;
0297   MonitorElement* meMatchQual_;
0298   MonitorElement* meDeltaZrealreal_;
0299   MonitorElement* meDeltaZfakefake_;
0300   MonitorElement* meDeltaZfakereal_;
0301   MonitorElement* meDeltaTrealreal_;
0302   MonitorElement* meDeltaTfakefake_;
0303   MonitorElement* meDeltaTfakereal_;
0304   MonitorElement* meRecoPosInSimCollection_;
0305   MonitorElement* meRecoPosInRecoOrigCollection_;
0306   MonitorElement* meSimPosInSimOrigCollection_;
0307   MonitorElement* meRecoPVPosSignal_;
0308   MonitorElement* meRecoPVPosSignalNotHighestPt_;
0309   MonitorElement* meRecoVtxVsLineDensity_;
0310   MonitorElement* meRecVerNumber_;
0311   MonitorElement* meRecPVZ_;
0312   MonitorElement* meRecPVT_;
0313   MonitorElement* meSimPVZ_;
0314 
0315   //some tests
0316   MonitorElement* meTrackResLowPTot_;
0317   MonitorElement* meTrackResHighPTot_;
0318   MonitorElement* meTrackPullLowPTot_;
0319   MonitorElement* meTrackPullHighPTot_;
0320 
0321   MonitorElement* meTrackResLowP_[3];
0322   MonitorElement* meTrackResHighP_[3];
0323   MonitorElement* meTrackPullLowP_[3];
0324   MonitorElement* meTrackPullHighP_[3];
0325 
0326   MonitorElement* meTrackResMassProtons_[3];
0327   MonitorElement* meTrackResMassTrueProtons_[3];
0328   MonitorElement* meTrackResMassPions_[3];
0329   MonitorElement* meTrackResMassTruePions_[3];
0330 
0331   MonitorElement* meBarrelPIDp_;
0332   MonitorElement* meEndcapPIDp_;
0333 
0334   MonitorElement* meBarrelNoPIDtype_;
0335   MonitorElement* meEndcapNoPIDtype_;
0336 
0337   MonitorElement* meBarrelTruePiNoPID_;
0338   MonitorElement* meBarrelTrueKNoPID_;
0339   MonitorElement* meBarrelTruePNoPID_;
0340   MonitorElement* meEndcapTruePiNoPID_;
0341   MonitorElement* meEndcapTrueKNoPID_;
0342   MonitorElement* meEndcapTruePNoPID_;
0343 
0344   MonitorElement* meBarrelTruePiAsPi_;
0345   MonitorElement* meBarrelTruePiAsK_;
0346   MonitorElement* meBarrelTruePiAsP_;
0347   MonitorElement* meEndcapTruePiAsPi_;
0348   MonitorElement* meEndcapTruePiAsK_;
0349   MonitorElement* meEndcapTruePiAsP_;
0350 
0351   MonitorElement* meBarrelTrueKAsPi_;
0352   MonitorElement* meBarrelTrueKAsK_;
0353   MonitorElement* meBarrelTrueKAsP_;
0354   MonitorElement* meEndcapTrueKAsPi_;
0355   MonitorElement* meEndcapTrueKAsK_;
0356   MonitorElement* meEndcapTrueKAsP_;
0357 
0358   MonitorElement* meBarrelTruePAsPi_;
0359   MonitorElement* meBarrelTruePAsK_;
0360   MonitorElement* meBarrelTruePAsP_;
0361   MonitorElement* meEndcapTruePAsPi_;
0362   MonitorElement* meEndcapTruePAsK_;
0363   MonitorElement* meEndcapTruePAsP_;
0364 };
0365 
0366 // constructors and destructor
0367 Primary4DVertexValidation::Primary4DVertexValidation(const edm::ParameterSet& iConfig)
0368     : folder_(iConfig.getParameter<std::string>("folder")),
0369       use_only_charged_tracks_(iConfig.getParameter<bool>("useOnlyChargedTracks")),
0370       debug_(iConfig.getUntrackedParameter<bool>("debug")),
0371       optionalPlots_(iConfig.getUntrackedParameter<bool>("optionalPlots")),
0372       minProbHeavy_(iConfig.getParameter<double>("minProbHeavy")),
0373       trackweightTh_(iConfig.getParameter<double>("trackweightTh")),
0374       mvaTh_(iConfig.getParameter<double>("mvaTh")),
0375       lineDensityPar_(iConfig.getParameter<std::vector<double>>("lineDensityPar")) {
0376   vecPileupSummaryInfoToken_ = consumes<std::vector<PileupSummaryInfo>>(edm::InputTag(std::string("addPileupInfo")));
0377   trackingParticleCollectionToken_ =
0378       consumes<TrackingParticleCollection>(iConfig.getParameter<edm::InputTag>("SimTag"));
0379   trackingVertexCollectionToken_ = consumes<TrackingVertexCollection>(iConfig.getParameter<edm::InputTag>("SimTag"));
0380   simToRecoAssociationToken_ =
0381       consumes<reco::SimToRecoCollection>(iConfig.getParameter<edm::InputTag>("TPtoRecoTrackAssoc"));
0382   recoToSimAssociationToken_ =
0383       consumes<reco::RecoToSimCollection>(iConfig.getParameter<edm::InputTag>("TPtoRecoTrackAssoc"));
0384   RecTrackToken_ = consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("mtdTracks"));
0385   RecBeamSpotToken_ = consumes<reco::BeamSpot>(iConfig.getParameter<edm::InputTag>("offlineBS"));
0386   Rec4DVerToken_ = consumes<edm::View<reco::Vertex>>(iConfig.getParameter<edm::InputTag>("offline4DPV"));
0387   trackAssocToken_ = consumes<edm::ValueMap<int>>(iConfig.getParameter<edm::InputTag>("trackAssocSrc"));
0388   pathLengthToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("pathLengthSrc"));
0389   momentumToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("momentumSrc"));
0390   timeToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("timeSrc"));
0391   t0PidToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("t0PID"));
0392   t0SafePidToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("t0SafePID"));
0393   sigmat0SafePidToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("sigmat0SafePID"));
0394   trackMVAQualToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("trackMVAQual"));
0395   tmtdToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("tmtd"));
0396   tofPiToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("tofPi"));
0397   tofKToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("tofK"));
0398   tofPToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("tofP"));
0399   probPiToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("probPi"));
0400   probKToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("probK"));
0401   probPToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("probP"));
0402 }
0403 
0404 Primary4DVertexValidation::~Primary4DVertexValidation() {}
0405 
0406 //
0407 // member functions
0408 //
0409 void Primary4DVertexValidation::bookHistograms(DQMStore::IBooker& ibook,
0410                                                edm::Run const& iRun,
0411                                                edm::EventSetup const& iSetup) {
0412   ibook.setCurrentFolder(folder_);
0413   // --- histograms booking
0414   meMVATrackEffPtTot_ = ibook.book1D("MVAEffPtTot", "Pt of tracks associated to LV; track pt [GeV] ", 110, 0., 11.);
0415   meMVATrackEffEtaTot_ = ibook.book1D("MVAEffEtaTot", "Pt of tracks associated to LV; track eta ", 66, 0., 3.3);
0416   meMVATrackMatchedEffPtTot_ =
0417       ibook.book1D("MVAMatchedEffPtTot", "Pt of tracks associated to LV matched to TP; track pt [GeV] ", 110, 0., 11.);
0418   meMVATrackMatchedEffPtMtd_ = ibook.book1D(
0419       "MVAMatchedEffPtMtd", "Pt of tracks associated to LV matched to TP with time; track pt [GeV] ", 110, 0., 11.);
0420   meMVATrackMatchedEffEtaTot_ =
0421       ibook.book1D("MVAMatchedEffEtaTot", "Pt of tracks associated to LV matched to TP; track eta ", 66, 0., 3.3);
0422   meMVATrackMatchedEffEtaMtd_ = ibook.book1D(
0423       "MVAMatchedEffEtaMtd", "Pt of tracks associated to LV matched to TP with time; track eta ", 66, 0., 3.3);
0424   meMVATrackResTot_ = ibook.book1D(
0425       "MVATrackRes", "t_{rec} - t_{sim} for tracks from LV MVA sel.; t_{rec} - t_{sim} [ns] ", 120, -0.15, 0.15);
0426   meTrackResTot_ = ibook.book1D("TrackRes", "t_{rec} - t_{sim} for tracks; t_{rec} - t_{sim} [ns] ", 120, -0.15, 0.15);
0427   meTrackRes_[0] = ibook.book1D(
0428       "TrackRes-LowMVA", "t_{rec} - t_{sim} for tracks with MVA < 0.5; t_{rec} - t_{sim} [ns] ", 100, -1., 1.);
0429   meTrackRes_[1] = ibook.book1D(
0430       "TrackRes-MediumMVA", "t_{rec} - t_{sim} for tracks with 0.5 < MVA < 0.8; t_{rec} - t_{sim} [ns] ", 100, -1., 1.);
0431   meTrackRes_[2] = ibook.book1D(
0432       "TrackRes-HighMVA", "t_{rec} - t_{sim} for tracks with MVA > 0.8; t_{rec} - t_{sim} [ns] ", 100, -1., 1.);
0433   if (optionalPlots_) {
0434     meTrackResMass_[0] = ibook.book1D(
0435         "TrackResMass-LowMVA", "t_{rec} - t_{est} for tracks with MVA < 0.5; t_{rec} - t_{est} [ns] ", 100, -1., 1.);
0436     meTrackResMass_[1] = ibook.book1D("TrackResMass-MediumMVA",
0437                                       "t_{rec} - t_{est} for tracks with 0.5 < MVA < 0.8; t_{rec} - t_{est} [ns] ",
0438                                       100,
0439                                       -1.,
0440                                       1.);
0441     meTrackResMass_[2] = ibook.book1D(
0442         "TrackResMass-HighMVA", "t_{rec} - t_{est} for tracks with MVA > 0.8; t_{rec} - t_{est} [ns] ", 100, -1., 1.);
0443     meTrackResMassTrue_[0] = ibook.book1D(
0444         "TrackResMassTrue-LowMVA", "t_{est} - t_{sim} for tracks with MVA < 0.5; t_{est} - t_{sim} [ns] ", 100, -1., 1.);
0445     meTrackResMassTrue_[1] = ibook.book1D("TrackResMassTrue-MediumMVA",
0446                                           "t_{est} - t_{sim} for tracks with 0.5 < MVA < 0.8; t_{est} - t_{sim} [ns] ",
0447                                           100,
0448                                           -1.,
0449                                           1.);
0450     meTrackResMassTrue_[2] = ibook.book1D("TrackResMassTrue-HighMVA",
0451                                           "t_{est} - t_{sim} for tracks with MVA > 0.8; t_{est} - t_{sim} [ns] ",
0452                                           100,
0453                                           -1.,
0454                                           1.);
0455   }
0456   meMVATrackPullTot_ =
0457       ibook.book1D("MVATrackPull", "Pull for tracks from LV MAV sel.; (t_{rec}-t_{sim})/#sigma_{t}", 50, -5., 5.);
0458   meTrackPullTot_ = ibook.book1D("TrackPull", "Pull for tracks; (t_{rec}-t_{sim})/#sigma_{t}", 100, -10., 10.);
0459   meTrackPull_[0] =
0460       ibook.book1D("TrackPull-LowMVA", "Pull for tracks with MVA < 0.5; (t_{rec}-t_{sim})/#sigma_{t}", 100, -10., 10.);
0461   meTrackPull_[1] = ibook.book1D(
0462       "TrackPull-MediumMVA", "Pull for tracks with 0.5 < MVA < 0.8; (t_{rec}-t_{sim})/#sigma_{t}", 100, -10., 10.);
0463   meTrackPull_[2] =
0464       ibook.book1D("TrackPull-HighMVA", "Pull for tracks with MVA > 0.8; (t_{rec}-t_{sim})/#sigma_{t}", 100, -10., 10.);
0465   meMVATrackZposResTot_ = ibook.book1D(
0466       "MVATrackZposResTot", "Z_{PCA} - Z_{sim} for tracks from LV MVA sel.;Z_{PCA} - Z_{sim} [cm] ", 50, -0.1, 0.1);
0467   meTrackZposResTot_ =
0468       ibook.book1D("TrackZposResTot", "Z_{PCA} - Z_{sim} for tracks;Z_{PCA} - Z_{sim} [cm] ", 50, -0.5, 0.5);
0469   meTrackZposRes_[0] = ibook.book1D(
0470       "TrackZposRes-LowMVA", "Z_{PCA} - Z_{sim} for tracks with MVA < 0.5;Z_{PCA} - Z_{sim} [cm] ", 50, -0.5, 0.5);
0471   meTrackZposRes_[1] = ibook.book1D("TrackZposRes-MediumMVA",
0472                                     "Z_{PCA} - Z_{sim} for tracks with 0.5 < MVA < 0.8 ;Z_{PCA} - Z_{sim} [cm] ",
0473                                     50,
0474                                     -0.5,
0475                                     0.5);
0476   meTrackZposRes_[2] = ibook.book1D(
0477       "TrackZposRes-HighMVA", "Z_{PCA} - Z_{sim} for tracks with MVA > 0.8 ;Z_{PCA} - Z_{sim} [cm] ", 50, -0.5, 0.5);
0478   meTrack3DposRes_[0] =
0479       ibook.book1D("Track3DposRes-LowMVA",
0480                    "3dPos_{PCA} - 3dPos_{sim} for tracks with MVA < 0.5 ;3dPos_{PCA} - 3dPos_{sim} [cm] ",
0481                    50,
0482                    -0.5,
0483                    0.5);
0484   meTrack3DposRes_[1] =
0485       ibook.book1D("Track3DposRes-MediumMVA",
0486                    "3dPos_{PCA} - 3dPos_{sim} for tracks with 0.5 < MVA < 0.8 ;3dPos_{PCA} - 3dPos_{sim} [cm] ",
0487                    50,
0488                    -0.5,
0489                    0.5);
0490   meTrack3DposRes_[2] =
0491       ibook.book1D("Track3DposRes-HighMVA",
0492                    "3dPos_{PCA} - 3dPos_{sim} for tracks with MVA > 0.8;3dPos_{PCA} - 3dPos_{sim} [cm] ",
0493                    50,
0494                    -0.5,
0495                    0.5);
0496   meTimeRes_ = ibook.book1D("TimeRes", "t_{rec} - t_{sim} ;t_{rec} - t_{sim} [ns] ", 40, -0.2, 0.2);
0497   meTimePull_ = ibook.book1D("TimePull", "Pull; t_{rec} - t_{sim}/#sigma_{t rec}", 100, -10., 10.);
0498   meTimeSignalRes_ =
0499       ibook.book1D("TimeSignalRes", "t_{rec} - t_{sim} for signal ;t_{rec} - t_{sim} [ns] ", 40, -0.2, 0.2);
0500   meTimeSignalPull_ =
0501       ibook.book1D("TimeSignalPull", "Pull for signal; t_{rec} - t_{sim}/#sigma_{t rec}", 100, -10., 10.);
0502   mePUvsRealV_ =
0503       ibook.bookProfile("PUvsReal", "#PU vertices vs #real matched vertices;#PU;#real ", 100, 0, 300, 100, 0, 200);
0504   mePUvsOtherFakeV_ = ibook.bookProfile(
0505       "PUvsOtherFake", "#PU vertices vs #other fake matched vertices;#PU;#other fake ", 100, 0, 300, 100, 0, 20);
0506   mePUvsSplitV_ =
0507       ibook.bookProfile("PUvsSplit", "#PU vertices vs #split matched vertices;#PU;#split ", 100, 0, 300, 100, 0, 20);
0508   meMatchQual_ = ibook.book1D("MatchQuality", "RECO-SIM vertex match quality; ", 8, 0, 8.);
0509   meDeltaZrealreal_ = ibook.book1D("DeltaZrealreal", "#Delta Z real-real; |#Delta Z (r-r)| [cm]", 100, 0, 0.5);
0510   meDeltaZfakefake_ = ibook.book1D("DeltaZfakefake", "#Delta Z fake-fake; |#Delta Z (f-f)| [cm]", 100, 0, 0.5);
0511   meDeltaZfakereal_ = ibook.book1D("DeltaZfakereal", "#Delta Z fake-real; |#Delta Z (f-r)| [cm]", 100, 0, 0.5);
0512   meDeltaTrealreal_ = ibook.book1D("DeltaTrealreal", "#Delta T real-real; |#Delta T (r-r)| [sigma]", 60, 0., 30.);
0513   meDeltaTfakefake_ = ibook.book1D("DeltaTfakefake", "#Delta T fake-fake; |#Delta T (f-f)| [sigma]", 60, 0., 30.);
0514   meDeltaTfakereal_ = ibook.book1D("DeltaTfakereal", "#Delta T fake-real; |#Delta T (f-r)| [sigma]", 60, 0., 30.);
0515   if (optionalPlots_) {
0516     meRecoPosInSimCollection_ = ibook.book1D(
0517         "RecoPosInSimCollection", "Sim signal vertex index associated to Reco signal vertex; Sim PV index", 200, 0, 200);
0518     meRecoPosInRecoOrigCollection_ =
0519         ibook.book1D("RecoPosInRecoOrigCollection", "Reco signal index in OrigCollection; Reco index", 200, 0, 200);
0520     meSimPosInSimOrigCollection_ =
0521         ibook.book1D("SimPosInSimOrigCollection", "Sim signal index in OrigCollection; Sim index", 200, 0, 200);
0522   }
0523   meRecoPVPosSignal_ =
0524       ibook.book1D("RecoPVPosSignal", "Position in reco collection of PV associated to sim signal", 200, 0, 200);
0525   meRecoPVPosSignalNotHighestPt_ =
0526       ibook.book1D("RecoPVPosSignalNotHighestPt",
0527                    "Position in reco collection of PV associated to sim signal not highest Pt",
0528                    200,
0529                    0,
0530                    200);
0531   meRecoVtxVsLineDensity_ =
0532       ibook.book1D("RecoVtxVsLineDensity", "#Reco vertices/mm/event; line density [#vtx/mm/event]", 160, 0., 4.);
0533   meRecVerNumber_ = ibook.book1D("RecVerNumber", "RECO Vertex Number: Number of vertices", 50, 0, 250);
0534   meRecPVZ_ = ibook.book1D("recPVZ", "Weighted #Rec vertices/mm", 400, -20., 20.);
0535   meRecPVT_ = ibook.book1D("recPVT", "#Rec vertices/10 ps", 200, -1., 1.);
0536   meSimPVZ_ = ibook.book1D("simPVZ", "Weighted #Sim vertices/mm", 400, -20., 20.);
0537 
0538   //some tests
0539   meTrackResLowPTot_ = ibook.book1D(
0540       "TrackResLowP", "t_{rec} - t_{sim} for tracks with p < 2 GeV; t_{rec} - t_{sim} [ns] ", 70, -0.15, 0.15);
0541   meTrackResLowP_[0] =
0542       ibook.book1D("TrackResLowP-LowMVA",
0543                    "t_{rec} - t_{sim} for tracks with MVA < 0.5 and p < 2 GeV; t_{rec} - t_{sim} [ns] ",
0544                    100,
0545                    -1.,
0546                    1.);
0547   meTrackResLowP_[1] =
0548       ibook.book1D("TrackResLowP-MediumMVA",
0549                    "t_{rec} - t_{sim} for tracks with 0.5 < MVA < 0.8 and p < 2 GeV; t_{rec} - t_{sim} [ns] ",
0550                    100,
0551                    -1.,
0552                    1.);
0553   meTrackResLowP_[2] =
0554       ibook.book1D("TrackResLowP-HighMVA",
0555                    "t_{rec} - t_{sim} for tracks with MVA > 0.8 and p < 2 GeV; t_{rec} - t_{sim} [ns] ",
0556                    100,
0557                    -1.,
0558                    1.);
0559   meTrackResHighPTot_ = ibook.book1D(
0560       "TrackResHighP", "t_{rec} - t_{sim} for tracks with p > 2 GeV; t_{rec} - t_{sim} [ns] ", 70, -0.15, 0.15);
0561   meTrackResHighP_[0] =
0562       ibook.book1D("TrackResHighP-LowMVA",
0563                    "t_{rec} - t_{sim} for tracks with MVA < 0.5 and p > 2 GeV; t_{rec} - t_{sim} [ns] ",
0564                    100,
0565                    -1.,
0566                    1.);
0567   meTrackResHighP_[1] =
0568       ibook.book1D("TrackResHighP-MediumMVA",
0569                    "t_{rec} - t_{sim} for tracks with 0.5 < MVA < 0.8 and p > 2 GeV; t_{rec} - t_{sim} [ns] ",
0570                    100,
0571                    -1.,
0572                    1.);
0573   meTrackResHighP_[2] =
0574       ibook.book1D("TrackResHighP-HighMVA",
0575                    "t_{rec} - t_{sim} for tracks with MVA > 0.8 and p > 2 GeV; t_{rec} - t_{sim} [ns] ",
0576                    100,
0577                    -1.,
0578                    1.);
0579   meTrackPullLowPTot_ =
0580       ibook.book1D("TrackPullLowP", "Pull for tracks with p < 2 GeV; (t_{rec}-t_{sim})/#sigma_{t}", 100, -10., 10.);
0581   meTrackPullLowP_[0] = ibook.book1D("TrackPullLowP-LowMVA",
0582                                      "Pull for tracks with MVA < 0.5 and p < 2 GeV; (t_{rec}-t_{sim})/#sigma_{t}",
0583                                      100,
0584                                      -10.,
0585                                      10.);
0586   meTrackPullLowP_[1] = ibook.book1D("TrackPullLowP-MediumMVA",
0587                                      "Pull for tracks with 0.5 < MVA < 0.8 and p < 2 GeV; (t_{rec}-t_{sim})/#sigma_{t}",
0588                                      100,
0589                                      -10.,
0590                                      10.);
0591   meTrackPullLowP_[2] = ibook.book1D("TrackPullLowP-HighMVA",
0592                                      "Pull for tracks with MVA > 0.8 and p < 2 GeV; (t_{rec}-t_{sim})/#sigma_{t}",
0593                                      100,
0594                                      -10.,
0595                                      10.);
0596   meTrackPullHighPTot_ =
0597       ibook.book1D("TrackPullHighP", "Pull for tracks with p > 2 GeV; (t_{rec}-t_{sim})/#sigma_{t}", 100, -10., 10.);
0598   meTrackPullHighP_[0] = ibook.book1D("TrackPullHighP-LowMVA",
0599                                       "Pull for tracks with MVA < 0.5 and p > 2 GeV; (t_{rec}-t_{sim})/#sigma_{t}",
0600                                       100,
0601                                       -10.,
0602                                       10.);
0603   meTrackPullHighP_[1] =
0604       ibook.book1D("TrackPullHighP-MediumMVA",
0605                    "Pull for tracks with 0.5 < MVA < 0.8 and p > 2 GeV; (t_{rec}-t_{sim})/#sigma_{t}",
0606                    100,
0607                    -10.,
0608                    10.);
0609   meTrackPullHighP_[2] = ibook.book1D("TrackPullHighP-HighMVA",
0610                                       "Pull for tracks with MVA > 0.8 and p > 2 GeV; (t_{rec}-t_{sim})/#sigma_{t}",
0611                                       100,
0612                                       -10.,
0613                                       10.);
0614   if (optionalPlots_) {
0615     meTrackResMassProtons_[0] =
0616         ibook.book1D("TrackResMass-Protons-LowMVA",
0617                      "t_{rec} - t_{est} for proton tracks with MVA < 0.5; t_{rec} - t_{est} [ns] ",
0618                      100,
0619                      -1.,
0620                      1.);
0621     meTrackResMassProtons_[1] =
0622         ibook.book1D("TrackResMass-Protons-MediumMVA",
0623                      "t_{rec} - t_{est} for proton tracks with 0.5 < MVA < 0.8; t_{rec} - t_{est} [ns] ",
0624                      100,
0625                      -1.,
0626                      1.);
0627     meTrackResMassProtons_[2] =
0628         ibook.book1D("TrackResMass-Protons-HighMVA",
0629                      "t_{rec} - t_{est} for proton tracks with MVA > 0.8; t_{rec} - t_{est} [ns] ",
0630                      100,
0631                      -1.,
0632                      1.);
0633     meTrackResMassTrueProtons_[0] =
0634         ibook.book1D("TrackResMassTrue-Protons-LowMVA",
0635                      "t_{est} - t_{sim} for proton tracks with MVA < 0.5; t_{est} - t_{sim} [ns] ",
0636                      100,
0637                      -1.,
0638                      1.);
0639     meTrackResMassTrueProtons_[1] =
0640         ibook.book1D("TrackResMassTrue-Protons-MediumMVA",
0641                      "t_{est} - t_{sim} for proton tracks with 0.5 < MVA < 0.8; t_{est} - t_{sim} [ns] ",
0642                      100,
0643                      -1.,
0644                      1.);
0645     meTrackResMassTrueProtons_[2] =
0646         ibook.book1D("TrackResMassTrue-Protons-HighMVA",
0647                      "t_{est} - t_{sim} for proton tracks with MVA > 0.8; t_{est} - t_{sim} [ns] ",
0648                      100,
0649                      -1.,
0650                      1.);
0651 
0652     meTrackResMassPions_[0] = ibook.book1D("TrackResMass-Pions-LowMVA",
0653                                            "t_{rec} - t_{est} for pion tracks with MVA < 0.5; t_{rec} - t_{est} [ns] ",
0654                                            100,
0655                                            -1.,
0656                                            1.);
0657     meTrackResMassPions_[1] =
0658         ibook.book1D("TrackResMass-Pions-MediumMVA",
0659                      "t_{rec} - t_{est} for pion tracks with 0.5 < MVA < 0.8; t_{rec} - t_{est} [ns] ",
0660                      100,
0661                      -1.,
0662                      1.);
0663     meTrackResMassPions_[2] = ibook.book1D("TrackResMass-Pions-HighMVA",
0664                                            "t_{rec} - t_{est} for pion tracks with MVA > 0.8; t_{rec} - t_{est} [ns] ",
0665                                            100,
0666                                            -1.,
0667                                            1.);
0668     meTrackResMassTruePions_[0] =
0669         ibook.book1D("TrackResMassTrue-Pions-LowMVA",
0670                      "t_{est} - t_{sim} for pion tracks with MVA < 0.5; t_{est} - t_{sim} [ns] ",
0671                      100,
0672                      -1.,
0673                      1.);
0674     meTrackResMassTruePions_[1] =
0675         ibook.book1D("TrackResMassTrue-Pions-MediumMVA",
0676                      "t_{est} - t_{sim} for pion tracks with 0.5 < MVA < 0.8; t_{est} - t_{sim} [ns] ",
0677                      100,
0678                      -1.,
0679                      1.);
0680     meTrackResMassTruePions_[2] =
0681         ibook.book1D("TrackResMassTrue-Pions-HighMVA",
0682                      "t_{est} - t_{sim} for pion tracks with MVA > 0.8; t_{est} - t_{sim} [ns] ",
0683                      100,
0684                      -1.,
0685                      1.);
0686   }
0687 
0688   meBarrelPIDp_ = ibook.book1D("BarrelPIDp", "PID track MTD momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
0689   meEndcapPIDp_ = ibook.book1D("EndcapPIDp", "PID track with MTD momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
0690 
0691   meBarrelNoPIDtype_ = ibook.book1D("BarrelNoPIDtype", "Barrel PID failure category", 4, 0., 4.);
0692   meEndcapNoPIDtype_ = ibook.book1D("EndcapNoPIDtype", "Endcap PID failure category", 4, 0., 4.);
0693 
0694   meBarrelNoPIDtype_ = ibook.book1D("BarrelNoPIDtype", "Barrel PID failure category", 4, 0., 4.);
0695   meEndcapNoPIDtype_ = ibook.book1D("EndcapNoPIDtype", "Endcap PID failure category", 4, 0., 4.);
0696 
0697   meBarrelTruePiNoPID_ =
0698       ibook.book1D("BarrelTruePiNoPID", "True pi NoPID momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
0699   meBarrelTrueKNoPID_ =
0700       ibook.book1D("BarrelTrueKNoPID", "True k NoPID momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
0701   meBarrelTruePNoPID_ =
0702       ibook.book1D("BarrelTruePNoPID", "True p NoPID momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
0703   meEndcapTruePiNoPID_ =
0704       ibook.book1D("EndcapTruePiNoPID", "True NoPIDpi momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
0705   meEndcapTrueKNoPID_ =
0706       ibook.book1D("EndcapTrueKNoPID", "True k NoPID momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
0707   meEndcapTruePNoPID_ =
0708       ibook.book1D("EndcapTruePNoPID", "True p NoPID momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
0709 
0710   meBarrelTruePiAsPi_ =
0711       ibook.book1D("BarrelTruePiAsPi", "True pi as pi momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
0712   meBarrelTruePiAsK_ =
0713       ibook.book1D("BarrelTruePiAsK", "True pi as k momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
0714   meBarrelTruePiAsP_ =
0715       ibook.book1D("BarrelTruePiAsP", "True pi as p momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
0716   meEndcapTruePiAsPi_ =
0717       ibook.book1D("EndcapTruePiAsPi", "True pi as pi momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
0718   meEndcapTruePiAsK_ =
0719       ibook.book1D("EndcapTruePiAsK", "True pi as k momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
0720   meEndcapTruePiAsP_ =
0721       ibook.book1D("EndcapTruePiAsP", "True pi as p momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
0722 
0723   meBarrelTrueKAsPi_ =
0724       ibook.book1D("BarrelTrueKAsPi", "True k as pi momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
0725   meBarrelTrueKAsK_ = ibook.book1D("BarrelTrueKAsK", "True k as k momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
0726   meBarrelTrueKAsP_ = ibook.book1D("BarrelTrueKAsP", "True k as p momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
0727   meEndcapTrueKAsPi_ =
0728       ibook.book1D("EndcapTrueKAsPi", "True k as pi momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
0729   meEndcapTrueKAsK_ = ibook.book1D("EndcapTrueKAsK", "True k as k momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
0730   meEndcapTrueKAsP_ = ibook.book1D("EndcapTrueKAsP", "True k as p momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
0731 
0732   meBarrelTruePAsPi_ =
0733       ibook.book1D("BarrelTruePAsPi", "True p as pi momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
0734   meBarrelTruePAsK_ = ibook.book1D("BarrelTruePAsK", "True p as k momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
0735   meBarrelTruePAsP_ = ibook.book1D("BarrelTruePAsP", "True p as p momentum spectrum, |eta| < 1.5;p [GeV]", 25, 0., 10.);
0736   meEndcapTruePAsPi_ =
0737       ibook.book1D("EndcapTruePAsPi", "True p as pi momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
0738   meEndcapTruePAsK_ = ibook.book1D("EndcapTruePAsK", "True p as k momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
0739   meEndcapTruePAsP_ = ibook.book1D("EndcapTruePAsP", "True p as p momentum spectrum, |eta| > 1.6;p [GeV]", 25, 0., 10.);
0740 }
0741 
0742 bool Primary4DVertexValidation::matchRecoTrack2SimSignal(const reco::TrackBaseRef& recoTrack) {
0743   auto found = r2s_->find(recoTrack);
0744 
0745   // reco track not matched to any TP
0746   if (found == r2s_->end())
0747     return false;
0748 
0749   //// reco track matched to some TP from signal vertex
0750   for (const auto& tp : found->val) {
0751     if (tp.first->eventId().bunchCrossing() == 0 && tp.first->eventId().event() == 0)
0752       return true;
0753   }
0754 
0755   // reco track not matched to any TP from signal vertex
0756   return false;
0757 }
0758 
0759 const edm::Ref<std::vector<TrackingParticle>>* Primary4DVertexValidation::getMatchedTP(
0760     const reco::TrackBaseRef& recoTrack, const TrackingVertexRef& vsim) {
0761   auto found = r2s_->find(recoTrack);
0762 
0763   // reco track not matched to any TP
0764   if (found == r2s_->end())
0765     return nullptr;
0766 
0767   //matched TP equal to any TP of sim vertex
0768   for (const auto& tp : found->val) {
0769     if (std::find_if(vsim->daughterTracks_begin(), vsim->daughterTracks_end(), [&](const TrackingParticleRef& vtp) {
0770           return tp.first == vtp;
0771         }) != vsim->daughterTracks_end())
0772       return &tp.first;
0773   }
0774 
0775   // reco track not matched to any TP from vertex
0776   return nullptr;
0777 }
0778 
0779 double Primary4DVertexValidation::timeFromTrueMass(double mass, double pathlength, double momentum, double time) {
0780   if (time > 0 && pathlength > 0 && mass > 0) {
0781     double gammasq = 1. + momentum * momentum / (mass * mass);
0782     double v = c_ * std::sqrt(1. - 1. / gammasq);  // cm / ns
0783     double t_est = time - (pathlength / v);
0784 
0785     return t_est;
0786   } else {
0787     return -1;
0788   }
0789 }
0790 
0791 bool Primary4DVertexValidation::select(const reco::Vertex& v, int level) {
0792   /* level
0793    0  !isFake  && ndof>4  (default)
0794    1  !isFake  && ndof>4 && prob > 0.01 
0795    2  !isFake  && ndof>4 && prob > 0.01 && ptmax2 > 0.4
0796    */
0797   if (v.isFake())
0798     return false;
0799   if ((level == 0) && (v.ndof() > selNdof_))
0800     return true;
0801   /*if ((level == 1) && (v.ndof() > selNdof_) && (vertex_pxy(v) > 0.01))
0802     return true;
0803   if ((level == 2) && (v.ndof() > selNdof_) && (vertex_pxy(v) > 0.01) && (vertex_ptmax2(v) > 0.4))
0804     return true;
0805   if ((level == 3) && (v.ndof() > selNdof_) && (vertex_ptmax2(v) < 0.4))
0806     return true;*/
0807   return false;
0808 }
0809 
0810 /* Extract information form TrackingParticles/TrackingVertex and fill
0811  * the helper class simPrimaryVertex with proper generation-level
0812  * information */
0813 std::vector<Primary4DVertexValidation::simPrimaryVertex> Primary4DVertexValidation::getSimPVs(
0814     const edm::Handle<TrackingVertexCollection>& tVC) {
0815   std::vector<Primary4DVertexValidation::simPrimaryVertex> simpv;
0816   int current_event = -1;
0817   int s = -1;
0818   for (TrackingVertexCollection::const_iterator v = tVC->begin(); v != tVC->end(); ++v) {
0819     //We keep only the first vertex from all the events at BX=0.
0820     if (v->eventId().bunchCrossing() != 0)
0821       continue;
0822     if (v->eventId().event() != current_event) {
0823       current_event = v->eventId().event();
0824     } else {
0825       continue;
0826     }
0827     s++;
0828     if (std::abs(v->position().z()) > 1000)
0829       continue;  // skip junk vertices
0830 
0831     // could be a new vertex, check  all primaries found so far to avoid multiple entries
0832     simPrimaryVertex sv(v->position().x(), v->position().y(), v->position().z(), v->position().t());
0833     sv.eventId = v->eventId();
0834     sv.sim_vertex = TrackingVertexRef(tVC, std::distance(tVC->begin(), v));
0835     sv.OriginalIndex = s;
0836 
0837     for (TrackingParticleRefVector::iterator iTrack = v->daughterTracks_begin(); iTrack != v->daughterTracks_end();
0838          ++iTrack) {
0839       assert((**iTrack).eventId().bunchCrossing() == 0);
0840     }
0841     simPrimaryVertex* vp = nullptr;  // will become non-NULL if a vertex is found and then point to it
0842     for (std::vector<simPrimaryVertex>::iterator v0 = simpv.begin(); v0 != simpv.end(); v0++) {
0843       if ((sv.eventId == v0->eventId) && (std::abs(sv.x - v0->x) < 1e-5) && (std::abs(sv.y - v0->y) < 1e-5) &&
0844           (std::abs(sv.z - v0->z) < 1e-5)) {
0845         vp = &(*v0);
0846         break;
0847       }
0848     }
0849     if (!vp) {
0850       // this is a new vertex, add it to the list of sim-vertices
0851       simpv.push_back(sv);
0852       vp = &simpv.back();
0853     }
0854 
0855     // Loop over daughter track(s) as Tracking Particles
0856     for (TrackingVertex::tp_iterator iTP = v->daughterTracks_begin(); iTP != v->daughterTracks_end(); ++iTP) {
0857       auto momentum = (*(*iTP)).momentum();
0858       const reco::Track* matched_best_reco_track = nullptr;
0859       double match_quality = -1;
0860       if (use_only_charged_tracks_ && (**iTP).charge() == 0)
0861         continue;
0862       if (s2r_->find(*iTP) != s2r_->end()) {
0863         matched_best_reco_track = (*s2r_)[*iTP][0].first.get();
0864         match_quality = (*s2r_)[*iTP][0].second;
0865       }
0866 
0867       vp->ptot.setPx(vp->ptot.x() + momentum.x());
0868       vp->ptot.setPy(vp->ptot.y() + momentum.y());
0869       vp->ptot.setPz(vp->ptot.z() + momentum.z());
0870       vp->ptot.setE(vp->ptot.e() + (**iTP).energy());
0871       vp->ptsq += ((**iTP).pt() * (**iTP).pt());
0872 
0873       if (matched_best_reco_track) {
0874         vp->num_matched_reco_tracks++;
0875         vp->average_match_quality += match_quality;
0876       }
0877     }  // End of for loop on daughters sim-particles
0878     if (vp->num_matched_reco_tracks)
0879       vp->average_match_quality /= static_cast<float>(vp->num_matched_reco_tracks);
0880     if (debug_) {
0881       edm::LogPrint("Primary4DVertexValidation")
0882           << "average number of associated tracks: " << vp->num_matched_reco_tracks / static_cast<float>(vp->nGenTrk)
0883           << " with average quality: " << vp->average_match_quality;
0884     }
0885   }  // End of for loop on tracking vertices
0886 
0887   // In case of no simulated vertices, break here
0888   if (simpv.empty())
0889     return simpv;
0890 
0891   // Now compute the closest distance in z between all simulated vertex
0892   // first initialize
0893   auto prev_z = simpv.back().z;
0894   for (simPrimaryVertex& vsim : simpv) {
0895     vsim.closest_vertex_distance_z = std::abs(vsim.z - prev_z);
0896     prev_z = vsim.z;
0897   }
0898   // then calculate
0899   for (std::vector<simPrimaryVertex>::iterator vsim = simpv.begin(); vsim != simpv.end(); vsim++) {
0900     std::vector<simPrimaryVertex>::iterator vsim2 = vsim;
0901     vsim2++;
0902     for (; vsim2 != simpv.end(); vsim2++) {
0903       double distance = std::abs(vsim->z - vsim2->z);
0904       // need both to be complete
0905       vsim->closest_vertex_distance_z = std::min(vsim->closest_vertex_distance_z, distance);
0906       vsim2->closest_vertex_distance_z = std::min(vsim2->closest_vertex_distance_z, distance);
0907     }
0908   }
0909   return simpv;
0910 }
0911 
0912 /* Extract information form recoVertex and fill the helper class
0913  * recoPrimaryVertex with proper reco-level information */
0914 std::vector<Primary4DVertexValidation::recoPrimaryVertex> Primary4DVertexValidation::getRecoPVs(
0915     const edm::Handle<edm::View<reco::Vertex>>& tVC) {
0916   std::vector<Primary4DVertexValidation::recoPrimaryVertex> recopv;
0917   int r = -1;
0918   for (auto v = tVC->begin(); v != tVC->end(); ++v) {
0919     r++;
0920     // Skip junk vertices
0921     if (std::abs(v->z()) > 1000)
0922       continue;
0923     if (v->isFake() || !v->isValid())
0924       continue;
0925 
0926     recoPrimaryVertex sv(v->position().x(), v->position().y(), v->position().z());
0927     sv.recVtx = &(*v);
0928     sv.recVtxRef = reco::VertexBaseRef(tVC, std::distance(tVC->begin(), v));
0929 
0930     sv.OriginalIndex = r;
0931     // this is a new vertex, add it to the list of reco-vertices
0932     recopv.push_back(sv);
0933     Primary4DVertexValidation::recoPrimaryVertex* vp = &recopv.back();
0934 
0935     // Loop over daughter track(s)
0936     for (auto iTrack = v->tracks_begin(); iTrack != v->tracks_end(); ++iTrack) {
0937       auto momentum = (*(*iTrack)).innerMomentum();
0938       if (momentum.mag2() == 0)
0939         momentum = (*(*iTrack)).momentum();
0940       vp->pt += std::sqrt(momentum.perp2());
0941       vp->ptsq += (momentum.perp2());
0942       vp->nRecoTrk++;
0943 
0944       auto matched = r2s_->find(*iTrack);
0945       if (matched != r2s_->end()) {
0946         vp->num_matched_sim_tracks++;
0947       }
0948 
0949     }  // End of for loop on daughters reconstructed tracks
0950   }    // End of for loop on tracking vertices
0951 
0952   // In case of no reco vertices, break here
0953   if (recopv.empty())
0954     return recopv;
0955 
0956   // Now compute the closest distance in z between all reconstructed vertex
0957   // first initialize
0958   auto prev_z = recopv.back().z;
0959   for (recoPrimaryVertex& vreco : recopv) {
0960     vreco.closest_vertex_distance_z = std::abs(vreco.z - prev_z);
0961     prev_z = vreco.z;
0962   }
0963   for (std::vector<recoPrimaryVertex>::iterator vreco = recopv.begin(); vreco != recopv.end(); vreco++) {
0964     std::vector<recoPrimaryVertex>::iterator vreco2 = vreco;
0965     vreco2++;
0966     for (; vreco2 != recopv.end(); vreco2++) {
0967       double distance = std::abs(vreco->z - vreco2->z);
0968       // need both to be complete
0969       vreco->closest_vertex_distance_z = std::min(vreco->closest_vertex_distance_z, distance);
0970       vreco2->closest_vertex_distance_z = std::min(vreco2->closest_vertex_distance_z, distance);
0971     }
0972   }
0973   return recopv;
0974 }
0975 
0976 // ------------ method called to produce the data  ------------
0977 void Primary4DVertexValidation::matchReco2Sim(std::vector<recoPrimaryVertex>& recopv,
0978                                               std::vector<simPrimaryVertex>& simpv,
0979                                               const edm::ValueMap<float>& sigmat0,
0980                                               const edm::ValueMap<float>& MVA,
0981                                               const edm::Handle<reco::BeamSpot>& BS) {
0982   for (auto vv : simpv) {
0983     vv.wnt.clear();
0984     vv.wos.clear();
0985   }
0986   for (auto rv : recopv) {
0987     rv.wnt.clear();
0988     rv.wos.clear();
0989   }
0990 
0991   for (unsigned int iv = 0; iv < recopv.size(); iv++) {
0992     const reco::Vertex* vertex = recopv.at(iv).recVtx;
0993 
0994     for (unsigned int iev = 0; iev < simpv.size(); iev++) {
0995       double wnt = 0;
0996       double wos = 0;
0997       double evwnt = 0;
0998       double evwos = 0;
0999       double evnt = 0;
1000 
1001       for (auto iTrack = vertex->tracks_begin(); iTrack != vertex->tracks_end(); ++iTrack) {
1002         double pt = (*iTrack)->pt();
1003 
1004         if (vertex->trackWeight(*iTrack) < trackweightTh_)
1005           continue;
1006         if (MVA[(*iTrack)] < mvaTh_)
1007           continue;
1008 
1009         auto tp_info = getMatchedTP(*iTrack, simpv.at(iev).sim_vertex);
1010         if (tp_info != nullptr) {
1011           double dz2_beam = pow((*BS).BeamWidthX() * cos((*iTrack)->phi()) / tan((*iTrack)->theta()), 2) +
1012                             pow((*BS).BeamWidthY() * sin((*iTrack)->phi()) / tan((*iTrack)->theta()), 2);
1013           double dz2 = pow((*iTrack)->dzError(), 2) + dz2_beam +
1014                        pow(0.0020, 2);  // added 20 um, some tracks have crazy small resolutions
1015           wos = vertex->trackWeight(*iTrack) / dz2;
1016           wnt = vertex->trackWeight(*iTrack) * std::min(pt, 1.0);
1017 
1018           if (sigmat0[(*iTrack)] > 0) {
1019             double sigmaZ = (*BS).sigmaZ();
1020             double sigmaT = sigmaZ / c_;  // c in cm/ns
1021             wos = wos / erf(sigmat0[(*iTrack)] / sigmaT);
1022           }
1023           simpv.at(iev).addTrack(iv, wos, wnt);
1024           recopv.at(iv).addTrack(iev, wos, wnt);
1025           evwos += wos;
1026           evwnt += wnt;
1027           evnt++;
1028         }
1029       }  //RecoTracks loop
1030 
1031       // require 2 tracks for a wos-match
1032       if ((evwos > 0) && (evwos > recopv.at(iv).maxwos) && (evnt > 1)) {
1033         recopv.at(iv).wosmatch = iev;
1034         recopv.at(iv).maxwos = evwos;
1035         recopv.at(iv).maxwosnt = evnt;
1036 
1037         simpv.at(iev).wos_dominated_recv.push_back(iv);
1038         simpv.at(iev).nwosmatch++;
1039       }
1040 
1041       // weighted track counting match, require at least one track
1042       if ((evnt > 0) && (evwnt > recopv.at(iv).maxwnt)) {
1043         recopv.at(iv).wntmatch = iev;
1044         recopv.at(iv).maxwnt = evwnt;
1045       }
1046     }  //TrackingVertex loop
1047 
1048   }  //RecoPrimaryVertex
1049 
1050   //after filling infos, goes for the sim-reco match
1051   for (auto& vrec : recopv) {
1052     vrec.sim = NOT_MATCHED;
1053     vrec.matchQuality = 0;
1054   }
1055   unsigned int iev = 0;
1056   for (auto& vv : simpv) {
1057     if (debug_) {
1058       edm::LogPrint("Primary4DVertexValidation") << "iev: " << iev;
1059       edm::LogPrint("Primary4DVertexValidation") << "wos_dominated_recv.size: " << vv.wos_dominated_recv.size();
1060     }
1061     for (unsigned int i = 0; i < vv.wos_dominated_recv.size(); i++) {
1062       auto recov = vv.wos_dominated_recv.at(i);
1063       if (debug_) {
1064         edm::LogPrint("Primary4DVertexValidation")
1065             << "index of reco vertex: " << recov << " that has a wos: " << vv.wos.at(recov) << " at position " << i;
1066       }
1067     }
1068     vv.rec = NOT_MATCHED;
1069     vv.matchQuality = 0;
1070     iev++;
1071   }
1072   //this tries a one-to-one match, taking simPV with highest wos if there are > 1 simPV candidates
1073   for (unsigned int rank = 1; rank < maxRank_; rank++) {
1074     for (unsigned int iev = 0; iev < simpv.size(); iev++) {  //loop on SimPV
1075       if (simpv.at(iev).rec != NOT_MATCHED)
1076         continue;
1077       if (simpv.at(iev).nwosmatch == 0)
1078         continue;
1079       if (simpv.at(iev).nwosmatch > rank)
1080         continue;
1081       unsigned int iv = NOT_MATCHED;
1082       for (unsigned int k = 0; k < simpv.at(iev).wos_dominated_recv.size(); k++) {
1083         unsigned int rec = simpv.at(iev).wos_dominated_recv.at(k);
1084         auto vrec = recopv.at(rec);
1085         if (vrec.sim != NOT_MATCHED)
1086           continue;  // already matched
1087         if (std::abs(simpv.at(iev).z - vrec.z) > zWosMatchMax_)
1088           continue;  // insanely far away
1089         if ((iv == NOT_MATCHED) || simpv.at(iev).wos.at(rec) > simpv.at(iev).wos.at(iv)) {
1090           iv = rec;
1091         }
1092       }
1093       if (iv !=
1094           NOT_MATCHED) {  //if the rec vertex has already been associated is possible that iv remains NOT_MATCHED at this point
1095         recopv.at(iv).sim = iev;
1096         simpv.at(iev).rec = iv;
1097         recopv.at(iv).matchQuality = rank;
1098         simpv.at(iev).matchQuality = rank;
1099       }
1100     }
1101   }
1102   //give vertices a chance that have a lot of overlap, but are still recognizably
1103   //caused by a specific simvertex (without being classified as dominating)
1104   //like a small peak sitting on the flank of a larger nearby peak
1105   unsigned int ntry = 0;
1106   while (ntry++ < maxTry_) {
1107     unsigned nmatch = 0;
1108     for (unsigned int iev = 0; iev < simpv.size(); iev++) {
1109       if ((simpv.at(iev).rec != NOT_MATCHED) || (simpv.at(iev).wos.empty()))
1110         continue;
1111       // find a rec vertex for the NOT_MATCHED sim vertex
1112       unsigned int rec = NOT_MATCHED;
1113       for (auto rv : simpv.at(iev).wos) {
1114         if ((rec == NOT_MATCHED) || (rv.second > simpv.at(iev).wos.at(rec))) {
1115           rec = rv.first;
1116         }
1117       }
1118 
1119       if (rec == NOT_MATCHED) {  //try with wnt match
1120         for (auto rv : simpv.at(iev).wnt) {
1121           if ((rec == NOT_MATCHED) || (rv.second > simpv.at(iev).wnt.at(rec))) {
1122             rec = rv.first;
1123           }
1124         }
1125       }
1126 
1127       if (rec == NOT_MATCHED)
1128         continue;
1129       if (recopv.at(rec).sim != NOT_MATCHED)
1130         continue;  // already gone
1131 
1132       // check if the recvertex can be  matched
1133       unsigned int rec2sim = NOT_MATCHED;
1134       for (auto sv : recopv.at(rec).wos) {
1135         if (simpv.at(sv.first).rec != NOT_MATCHED)
1136           continue;  // already used
1137         if ((rec2sim == NOT_MATCHED) || (sv.second > recopv.at(rec).wos.at(rec2sim))) {
1138           rec2sim = sv.first;
1139         }
1140       }
1141       if (iev == rec2sim) {
1142         // do the match and assign lowest quality (i.e. max rank)
1143         recopv.at(rec).sim = iev;
1144         recopv.at(rec).matchQuality = maxRank_;
1145         simpv.at(iev).rec = rec;
1146         simpv.at(iev).matchQuality = maxRank_;
1147         nmatch++;
1148       }
1149     }  //sim loop
1150     if (nmatch == 0) {
1151       break;
1152     }
1153   }  // ntry
1154 }
1155 
1156 void Primary4DVertexValidation::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
1157   using edm::Handle;
1158   using edm::View;
1159   using std::cout;
1160   using std::endl;
1161   using std::vector;
1162   using namespace reco;
1163 
1164   std::vector<float> pileUpInfo_z;
1165 
1166   // get the pileup information
1167   edm::Handle<std::vector<PileupSummaryInfo>> puinfoH;
1168   if (iEvent.getByToken(vecPileupSummaryInfoToken_, puinfoH)) {
1169     for (auto const& pu_info : *puinfoH.product()) {
1170       if (pu_info.getBunchCrossing() == 0) {
1171         pileUpInfo_z = pu_info.getPU_zpositions();
1172         break;
1173       }
1174     }
1175   }
1176 
1177   edm::Handle<TrackingParticleCollection> TPCollectionH;
1178   iEvent.getByToken(trackingParticleCollectionToken_, TPCollectionH);
1179   if (!TPCollectionH.isValid())
1180     edm::LogWarning("Primary4DVertexValidation") << "TPCollectionH is not valid";
1181 
1182   edm::Handle<TrackingVertexCollection> TVCollectionH;
1183   iEvent.getByToken(trackingVertexCollectionToken_, TVCollectionH);
1184   if (!TVCollectionH.isValid())
1185     edm::LogWarning("Primary4DVertexValidation") << "TVCollectionH is not valid";
1186 
1187   edm::Handle<reco::SimToRecoCollection> simToRecoH;
1188   iEvent.getByToken(simToRecoAssociationToken_, simToRecoH);
1189   if (simToRecoH.isValid())
1190     s2r_ = simToRecoH.product();
1191   else
1192     edm::LogWarning("Primary4DVertexValidation") << "simToRecoH is not valid";
1193 
1194   edm::Handle<reco::RecoToSimCollection> recoToSimH;
1195   iEvent.getByToken(recoToSimAssociationToken_, recoToSimH);
1196   if (recoToSimH.isValid())
1197     r2s_ = recoToSimH.product();
1198   else
1199     edm::LogWarning("Primary4DVertexValidation") << "recoToSimH is not valid";
1200 
1201   edm::Handle<reco::BeamSpot> BeamSpotH;
1202   iEvent.getByToken(RecBeamSpotToken_, BeamSpotH);
1203   if (!BeamSpotH.isValid())
1204     edm::LogWarning("Primary4DVertexValidation") << "BeamSpotH is not valid";
1205 
1206   std::vector<simPrimaryVertex> simpv;  // a list of simulated primary MC vertices
1207   simpv = getSimPVs(TVCollectionH);
1208   //this bool check if first vertex in that with highest pT
1209   bool signal_is_highest_pt =
1210       std::max_element(simpv.begin(), simpv.end(), [](const simPrimaryVertex& lhs, const simPrimaryVertex& rhs) {
1211         return lhs.ptsq < rhs.ptsq;
1212       }) == simpv.begin();
1213 
1214   std::vector<recoPrimaryVertex> recopv;  // a list of reconstructed primary MC vertices
1215   edm::Handle<edm::View<reco::Vertex>> recVtxs;
1216   iEvent.getByToken(Rec4DVerToken_, recVtxs);
1217   if (!recVtxs.isValid())
1218     edm::LogWarning("Primary4DVertexValidation") << "recVtxs is not valid";
1219   recopv = getRecoPVs(recVtxs);
1220 
1221   const auto& trackAssoc = iEvent.get(trackAssocToken_);
1222   const auto& pathLength = iEvent.get(pathLengthToken_);
1223   const auto& momentum = iEvent.get(momentumToken_);
1224   const auto& time = iEvent.get(timeToken_);
1225   const auto& t0Pid = iEvent.get(t0PidToken_);
1226   const auto& t0Safe = iEvent.get(t0SafePidToken_);
1227   const auto& sigmat0Safe = iEvent.get(sigmat0SafePidToken_);
1228   const auto& mtdQualMVA = iEvent.get(trackMVAQualToken_);
1229   const auto& tMtd = iEvent.get(tmtdToken_);
1230   const auto& tofPi = iEvent.get(tofPiToken_);
1231   const auto& tofK = iEvent.get(tofKToken_);
1232   const auto& tofP = iEvent.get(tofPToken_);
1233   const auto& probPi = iEvent.get(probPiToken_);
1234   const auto& probK = iEvent.get(probKToken_);
1235   const auto& probP = iEvent.get(probPToken_);
1236 
1237   //I have simPV and recoPV collections
1238   matchReco2Sim(recopv, simpv, sigmat0Safe, mtdQualMVA, BeamSpotH);
1239 
1240   //Loop on tracks
1241   for (unsigned int iv = 0; iv < recopv.size(); iv++) {
1242     const reco::Vertex* vertex = recopv.at(iv).recVtx;
1243 
1244     for (unsigned int iev = 0; iev < simpv.size(); iev++) {
1245       auto vsim = simpv.at(iev).sim_vertex;
1246 
1247       bool selectedVtxMatching = recopv.at(iv).sim == iev && simpv.at(iev).rec == iv &&
1248                                  simpv.at(iev).eventId.bunchCrossing() == 0 && simpv.at(iev).eventId.event() == 0 &&
1249                                  recopv.at(iv).OriginalIndex == 0;
1250       if (selectedVtxMatching && !recopv.at(iv).is_signal()) {
1251         edm::LogWarning("Primary4DVertexValidation")
1252             << "Reco vtx leading match inconsistent: BX/ID " << simpv.at(iev).eventId.bunchCrossing() << " "
1253             << simpv.at(iev).eventId.event();
1254       }
1255       double vzsim = simpv.at(iev).z;
1256       double vtsim = simpv.at(iev).t * simUnit_;
1257 
1258       for (auto iTrack = vertex->tracks_begin(); iTrack != vertex->tracks_end(); ++iTrack) {
1259         if (trackAssoc[*iTrack] == -1) {
1260           LogTrace("mtdTracks") << "Extended track not associated";
1261           continue;
1262         }
1263 
1264         if (vertex->trackWeight(*iTrack) < trackweightTh_)
1265           continue;
1266 
1267         bool noCrack = std::abs((*iTrack)->eta()) < trackMaxBtlEta_ || std::abs((*iTrack)->eta()) > trackMinEtlEta_;
1268 
1269         bool selectRecoTrk = mvaRecSel(**iTrack, *vertex, t0Safe[*iTrack], sigmat0Safe[*iTrack]);
1270         if (selectedVtxMatching && selectRecoTrk) {
1271           if (noCrack) {
1272             meMVATrackEffPtTot_->Fill((*iTrack)->pt());
1273           }
1274           meMVATrackEffEtaTot_->Fill(std::abs((*iTrack)->eta()));
1275         }
1276 
1277         auto tp_info = getMatchedTP(*iTrack, vsim);
1278         if (tp_info != nullptr) {
1279           double mass = (*tp_info)->mass();
1280           double tsim = (*tp_info)->parentVertex()->position().t() * simUnit_;
1281           double tEst = timeFromTrueMass(mass, pathLength[*iTrack], momentum[*iTrack], time[*iTrack]);
1282 
1283           double xsim = (*tp_info)->parentVertex()->position().x();
1284           double ysim = (*tp_info)->parentVertex()->position().y();
1285           double zsim = (*tp_info)->parentVertex()->position().z();
1286           double xPCA = (*iTrack)->vx();
1287           double yPCA = (*iTrack)->vy();
1288           double zPCA = (*iTrack)->vz();
1289 
1290           double dZ = zPCA - zsim;
1291           double d3D = std::sqrt((xPCA - xsim) * (xPCA - xsim) + (yPCA - ysim) * (yPCA - ysim) + dZ * dZ);
1292           // orient d3D according to the projection of RECO - SIM onto simulated momentum
1293           if ((xPCA - xsim) * ((*tp_info)->px()) + (yPCA - ysim) * ((*tp_info)->py()) + dZ * ((*tp_info)->pz()) < 0.) {
1294             d3D = -d3D;
1295           }
1296 
1297           bool selectTP = mvaTPSel(**tp_info);
1298 
1299           if (selectedVtxMatching && selectRecoTrk && selectTP) {
1300             meMVATrackZposResTot_->Fill((*iTrack)->vz() - vzsim);
1301             if (noCrack) {
1302               meMVATrackMatchedEffPtTot_->Fill((*iTrack)->pt());
1303             }
1304             meMVATrackMatchedEffEtaTot_->Fill(std::abs((*iTrack)->eta()));
1305           }
1306 
1307           if (sigmat0Safe[*iTrack] == -1)
1308             continue;
1309 
1310           if (selectedVtxMatching && selectRecoTrk && selectTP) {
1311             meMVATrackResTot_->Fill(t0Safe[*iTrack] - vtsim);
1312             meMVATrackPullTot_->Fill((t0Safe[*iTrack] - vtsim) / sigmat0Safe[*iTrack]);
1313             if (noCrack) {
1314               meMVATrackMatchedEffPtMtd_->Fill((*iTrack)->pt());
1315             }
1316             meMVATrackMatchedEffEtaMtd_->Fill(std::abs((*iTrack)->eta()));
1317 
1318             unsigned int noPIDtype = 0;
1319             if (probPi[*iTrack] == -1) {
1320               noPIDtype = 1;
1321             } else if (isnan(probPi[*iTrack])) {
1322               noPIDtype = 2;
1323             } else if (probPi[*iTrack] == 1 && probK[*iTrack] == 0 && probP[*iTrack] == 0) {
1324               noPIDtype = 3;
1325             }
1326             bool noPID = noPIDtype > 0;
1327             bool isPi = !noPID && 1. - probPi[*iTrack] < minProbHeavy_;
1328             bool isK = !noPID && !isPi && probK[*iTrack] > probP[*iTrack];
1329             bool isP = !noPID && !isPi && !isK;
1330 
1331             if ((isPi && std::abs(tMtd[*iTrack] - tofPi[*iTrack] - t0Pid[*iTrack]) > tol_) ||
1332                 (isK && std::abs(tMtd[*iTrack] - tofK[*iTrack] - t0Pid[*iTrack]) > tol_) ||
1333                 (isP && std::abs(tMtd[*iTrack] - tofP[*iTrack] - t0Pid[*iTrack]) > tol_)) {
1334               edm::LogWarning("Primary4DVertexValidation")
1335                   << "No match between mass hyp. and time: " << std::abs((*tp_info)->pdgId()) << " mass hyp pi/k/p "
1336                   << isPi << " " << isK << " " << isP << " t0/t0safe " << t0Pid[*iTrack] << " " << t0Safe[*iTrack]
1337                   << " tMtd - tof pi/K/p " << tMtd[*iTrack] - tofPi[*iTrack] << " " << tMtd[*iTrack] - tofK[*iTrack]
1338                   << " " << tMtd[*iTrack] - tofP[*iTrack] << " Prob pi/K/p " << probPi[*iTrack] << " " << probK[*iTrack]
1339                   << " " << probP[*iTrack];
1340             }
1341 
1342             if (std::abs((*iTrack)->eta()) < trackMaxBtlEta_) {
1343               meBarrelPIDp_->Fill((*iTrack)->p());
1344               meBarrelNoPIDtype_->Fill(noPIDtype + 0.5);
1345               if (std::abs((*tp_info)->pdgId()) == 211) {
1346                 if (noPID) {
1347                   meBarrelTruePiNoPID_->Fill((*iTrack)->p());
1348                 } else if (isPi) {
1349                   meBarrelTruePiAsPi_->Fill((*iTrack)->p());
1350                 } else if (isK) {
1351                   meBarrelTruePiAsK_->Fill((*iTrack)->p());
1352                 } else if (isP) {
1353                   meBarrelTruePiAsP_->Fill((*iTrack)->p());
1354                 } else {
1355                   edm::LogWarning("Primary4DVertexValidation")
1356                       << "No PID class: " << std::abs((*tp_info)->pdgId()) << " t0/t0safe " << t0Pid[*iTrack] << " "
1357                       << t0Safe[*iTrack] << " Prob pi/K/p " << probPi[*iTrack] << " " << probK[*iTrack] << " "
1358                       << probP[*iTrack];
1359                 }
1360               } else if (std::abs((*tp_info)->pdgId()) == 321) {
1361                 if (noPID) {
1362                   meBarrelTrueKNoPID_->Fill((*iTrack)->p());
1363                 } else if (isPi) {
1364                   meBarrelTrueKAsPi_->Fill((*iTrack)->p());
1365                 } else if (isK) {
1366                   meBarrelTrueKAsK_->Fill((*iTrack)->p());
1367                 } else if (isP) {
1368                   meBarrelTrueKAsP_->Fill((*iTrack)->p());
1369                 } else {
1370                   edm::LogWarning("Primary4DVertexValidation")
1371                       << "No PID class: " << std::abs((*tp_info)->pdgId()) << " t0/t0safe " << t0Pid[*iTrack] << " "
1372                       << t0Safe[*iTrack] << " Prob pi/K/p " << probPi[*iTrack] << " " << probK[*iTrack] << " "
1373                       << probP[*iTrack];
1374                 }
1375               } else if (std::abs((*tp_info)->pdgId()) == 2212) {
1376                 if (noPID) {
1377                   meBarrelTruePNoPID_->Fill((*iTrack)->p());
1378                 } else if (isPi) {
1379                   meBarrelTruePAsPi_->Fill((*iTrack)->p());
1380                 } else if (isK) {
1381                   meBarrelTruePAsK_->Fill((*iTrack)->p());
1382                 } else if (isP) {
1383                   meBarrelTruePAsP_->Fill((*iTrack)->p());
1384                 } else {
1385                   edm::LogWarning("Primary4DVertexValidation")
1386                       << "No PID class: " << std::abs((*tp_info)->pdgId()) << " t0/t0safe " << t0Pid[*iTrack] << " "
1387                       << t0Safe[*iTrack] << " Prob pi/K/p " << probPi[*iTrack] << " " << probK[*iTrack] << " "
1388                       << probP[*iTrack];
1389                 }
1390               }
1391             } else if (std::abs((*iTrack)->eta()) > trackMinEtlEta_ && std::abs((*iTrack)->eta()) < trackMaxEtlEta_) {
1392               meEndcapPIDp_->Fill((*iTrack)->p());
1393               meEndcapNoPIDtype_->Fill(noPIDtype + 0.5);
1394               if (std::abs((*tp_info)->pdgId()) == 211) {
1395                 if (noPID) {
1396                   meEndcapTruePiNoPID_->Fill((*iTrack)->p());
1397                 } else if (isPi) {
1398                   meEndcapTruePiAsPi_->Fill((*iTrack)->p());
1399                 } else if (isK) {
1400                   meEndcapTruePiAsK_->Fill((*iTrack)->p());
1401                 } else if (isP) {
1402                   meEndcapTruePiAsP_->Fill((*iTrack)->p());
1403                 } else {
1404                   edm::LogWarning("Primary4DVertexValidation")
1405                       << "No PID class: " << std::abs((*tp_info)->pdgId()) << " t0/t0safe " << t0Pid[*iTrack] << " "
1406                       << t0Safe[*iTrack] << " Prob pi/K/p " << probPi[*iTrack] << " " << probK[*iTrack] << " "
1407                       << probP[*iTrack];
1408                 }
1409               } else if (std::abs((*tp_info)->pdgId()) == 321) {
1410                 if (noPID) {
1411                   meEndcapTrueKNoPID_->Fill((*iTrack)->p());
1412                 } else if (isPi) {
1413                   meEndcapTrueKAsPi_->Fill((*iTrack)->p());
1414                 } else if (isK) {
1415                   meEndcapTrueKAsK_->Fill((*iTrack)->p());
1416                 } else if (isP) {
1417                   meEndcapTrueKAsP_->Fill((*iTrack)->p());
1418                 } else {
1419                   edm::LogWarning("Primary4DVertexValidation")
1420                       << "No PID class: " << std::abs((*tp_info)->pdgId()) << " t0/t0safe " << t0Pid[*iTrack] << " "
1421                       << t0Safe[*iTrack] << " Prob pi/K/p " << probPi[*iTrack] << " " << probK[*iTrack] << " "
1422                       << probP[*iTrack];
1423                 }
1424               } else if (std::abs((*tp_info)->pdgId()) == 2212) {
1425                 if (noPID) {
1426                   meEndcapTruePNoPID_->Fill((*iTrack)->p());
1427                 } else if (isPi) {
1428                   meEndcapTruePAsPi_->Fill((*iTrack)->p());
1429                 } else if (isK) {
1430                   meEndcapTruePAsK_->Fill((*iTrack)->p());
1431                 } else if (isP) {
1432                   meEndcapTruePAsP_->Fill((*iTrack)->p());
1433                 } else {
1434                   edm::LogWarning("Primary4DVertexValidation")
1435                       << "No PID class: " << std::abs((*tp_info)->pdgId()) << " t0/t0safe " << t0Pid[*iTrack] << " "
1436                       << t0Safe[*iTrack] << " Prob pi/K/p " << probPi[*iTrack] << " " << probK[*iTrack] << " "
1437                       << probP[*iTrack];
1438                 }
1439               }
1440             }
1441           }
1442           meTrackResTot_->Fill(t0Safe[*iTrack] - tsim);
1443           meTrackPullTot_->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
1444           meTrackZposResTot_->Fill(dZ);
1445           if ((*iTrack)->p() <= 2) {
1446             meTrackResLowPTot_->Fill(t0Safe[*iTrack] - tsim);
1447             meTrackPullLowPTot_->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
1448           } else {
1449             meTrackResHighPTot_->Fill(t0Safe[*iTrack] - tsim);
1450             meTrackPullHighPTot_->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
1451           }
1452 
1453           if (mtdQualMVA[(*iTrack)] < mvaL_) {
1454             meTrackZposRes_[0]->Fill(dZ);
1455             meTrack3DposRes_[0]->Fill(d3D);
1456             meTrackRes_[0]->Fill(t0Safe[*iTrack] - tsim);
1457             meTrackPull_[0]->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
1458 
1459             if (optionalPlots_) {
1460               meTrackResMass_[0]->Fill(t0Safe[*iTrack] - tEst);
1461               meTrackResMassTrue_[0]->Fill(tEst - tsim);
1462             }
1463 
1464             if ((*iTrack)->p() <= 2) {
1465               meTrackResLowP_[0]->Fill(t0Safe[*iTrack] - tsim);
1466               meTrackPullLowP_[0]->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
1467             } else if ((*iTrack)->p() > 2) {
1468               meTrackResHighP_[0]->Fill(t0Safe[*iTrack] - tsim);
1469               meTrackPullHighP_[0]->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
1470             }
1471 
1472             if (optionalPlots_) {
1473               if (std::abs((*tp_info)->pdgId()) == 2212) {
1474                 meTrackResMassProtons_[0]->Fill(t0Safe[*iTrack] - tEst);
1475                 meTrackResMassTrueProtons_[0]->Fill(tEst - tsim);
1476               } else if (std::abs((*tp_info)->pdgId()) == 211) {
1477                 meTrackResMassPions_[0]->Fill(t0Safe[*iTrack] - tEst);
1478                 meTrackResMassTruePions_[0]->Fill(tEst - tsim);
1479               }
1480             }
1481 
1482           } else if (mtdQualMVA[(*iTrack)] > mvaL_ && mtdQualMVA[(*iTrack)] < mvaH_) {
1483             meTrackZposRes_[1]->Fill(dZ);
1484             meTrack3DposRes_[1]->Fill(d3D);
1485             meTrackRes_[1]->Fill(t0Safe[*iTrack] - tsim);
1486             meTrackPull_[1]->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
1487 
1488             if (optionalPlots_) {
1489               meTrackResMass_[1]->Fill(t0Safe[*iTrack] - tEst);
1490               meTrackResMassTrue_[1]->Fill(tEst - tsim);
1491             }
1492 
1493             if ((*iTrack)->p() <= 2) {
1494               meTrackResLowP_[1]->Fill(t0Safe[*iTrack] - tsim);
1495               meTrackPullLowP_[1]->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
1496             } else if ((*iTrack)->p() > 2) {
1497               meTrackResHighP_[1]->Fill(t0Safe[*iTrack] - tsim);
1498               meTrackPullHighP_[1]->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
1499             }
1500 
1501             if (optionalPlots_) {
1502               if (std::abs((*tp_info)->pdgId()) == 2212) {
1503                 meTrackResMassProtons_[1]->Fill(t0Safe[*iTrack] - tEst);
1504                 meTrackResMassTrueProtons_[1]->Fill(tEst - tsim);
1505               } else if (std::abs((*tp_info)->pdgId()) == 211) {
1506                 meTrackResMassPions_[1]->Fill(t0Safe[*iTrack] - tEst);
1507                 meTrackResMassTruePions_[1]->Fill(tEst - tsim);
1508               }
1509             }
1510 
1511           } else if (mtdQualMVA[(*iTrack)] > mvaH_) {
1512             meTrackZposRes_[2]->Fill(dZ);
1513             meTrack3DposRes_[2]->Fill(d3D);
1514             meTrackRes_[2]->Fill(t0Safe[*iTrack] - tsim);
1515             meTrackPull_[2]->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
1516 
1517             if (optionalPlots_) {
1518               meTrackResMass_[2]->Fill(t0Safe[*iTrack] - tEst);
1519               meTrackResMassTrue_[2]->Fill(tEst - tsim);
1520             }
1521 
1522             if ((*iTrack)->p() <= 2) {
1523               meTrackResLowP_[2]->Fill(t0Safe[*iTrack] - tsim);
1524               meTrackPullLowP_[2]->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
1525             } else if ((*iTrack)->p() > 2) {
1526               meTrackResHighP_[2]->Fill(t0Safe[*iTrack] - tsim);
1527               meTrackPullHighP_[2]->Fill((t0Safe[*iTrack] - tsim) / sigmat0Safe[*iTrack]);
1528             }
1529 
1530             if (optionalPlots_) {
1531               if (std::abs((*tp_info)->pdgId()) == 2212) {
1532                 meTrackResMassProtons_[2]->Fill(t0Safe[*iTrack] - tEst);
1533                 meTrackResMassTrueProtons_[2]->Fill(tEst - tsim);
1534               } else if (std::abs((*tp_info)->pdgId()) == 211) {
1535                 meTrackResMassPions_[2]->Fill(t0Safe[*iTrack] - tEst);
1536                 meTrackResMassTruePions_[2]->Fill(tEst - tsim);
1537               }
1538             }
1539           }
1540         }  //if tp_info != nullptr
1541       }
1542     }
1543   }
1544 
1545   int real = 0;
1546   int fake = 0;
1547   int other_fake = 0;
1548   int split = 0;
1549 
1550   auto puLineDensity = [&](double z) {
1551     // gaussian parameterization of line density vs z, z in cm, parameters in mm
1552     double argl = (z * 10. - lineDensityPar_[1]) / lineDensityPar_[2];
1553     return lineDensityPar_[0] * exp(-0.5 * argl * argl);
1554   };
1555 
1556   meRecVerNumber_->Fill(recopv.size());
1557   for (unsigned int ir = 0; ir < recopv.size(); ir++) {
1558     meRecoVtxVsLineDensity_->Fill(puLineDensity(recopv.at(ir).z));
1559     meRecPVZ_->Fill(recopv.at(ir).z, 1. / puLineDensity(recopv.at(ir).z));
1560     if (recopv.at(ir).recVtx->tError() > 0.) {
1561       meRecPVT_->Fill(recopv.at(ir).recVtx->t());
1562     }
1563     if (debug_) {
1564       edm::LogPrint("Primary4DVertexValidation") << "************* IR: " << ir;
1565       edm::LogPrint("Primary4DVertexValidation")
1566           << "z: " << recopv.at(ir).z << " corresponding to line density: " << puLineDensity(recopv.at(ir).z);
1567       edm::LogPrint("Primary4DVertexValidation") << "is_real: " << recopv.at(ir).is_real();
1568       edm::LogPrint("Primary4DVertexValidation") << "is_fake: " << recopv.at(ir).is_fake();
1569       edm::LogPrint("Primary4DVertexValidation") << "is_signal: " << recopv.at(ir).is_signal();
1570       edm::LogPrint("Primary4DVertexValidation") << "split_from: " << recopv.at(ir).split_from();
1571       edm::LogPrint("Primary4DVertexValidation") << "other fake: " << recopv.at(ir).other_fake();
1572     }
1573     if (recopv.at(ir).is_real())
1574       real++;
1575     if (recopv.at(ir).is_fake())
1576       fake++;
1577     if (recopv.at(ir).other_fake())
1578       other_fake++;
1579     if (recopv.at(ir).split_from() != -1) {
1580       split++;
1581     }
1582   }
1583 
1584   if (debug_) {
1585     edm::LogPrint("Primary4DVertexValidation") << "is_real: " << real;
1586     edm::LogPrint("Primary4DVertexValidation") << "is_fake: " << fake;
1587     edm::LogPrint("Primary4DVertexValidation") << "split_from: " << split;
1588     edm::LogPrint("Primary4DVertexValidation") << "other fake: " << other_fake;
1589   }
1590 
1591   //fill vertices histograms here in a new loop
1592   for (unsigned int is = 0; is < simpv.size(); is++) {
1593     // protect against particle guns with very displaced vertices
1594     if (std::isinf(1. / puLineDensity(simpv.at(is).z))) {
1595       continue;
1596     }
1597     meSimPVZ_->Fill(simpv.at(is).z, 1. / puLineDensity(simpv.at(is).z));
1598     if (is == 0 && optionalPlots_) {
1599       meSimPosInSimOrigCollection_->Fill(simpv.at(is).OriginalIndex);
1600     }
1601 
1602     if (simpv.at(is).rec == NOT_MATCHED) {
1603       if (debug_) {
1604         edm::LogPrint("Primary4DVertexValidation") << "sim vertex: " << is << " is not matched with any reco";
1605       }
1606       continue;
1607     }
1608 
1609     for (unsigned int ir = 0; ir < recopv.size(); ir++) {
1610       if (recopv.at(ir).sim == is && simpv.at(is).rec == ir) {
1611         meTimeRes_->Fill(recopv.at(ir).recVtx->t() - simpv.at(is).t * simUnit_);
1612         meTimePull_->Fill((recopv.at(ir).recVtx->t() - simpv.at(is).t * simUnit_) / recopv.at(ir).recVtx->tError());
1613         mePUvsRealV_->Fill(simpv.size(), real);
1614         mePUvsOtherFakeV_->Fill(simpv.size(), other_fake);
1615         mePUvsSplitV_->Fill(simpv.size(), split);
1616         meMatchQual_->Fill(recopv.at(ir).matchQuality - 0.5);
1617         if (ir == 0) {  //signal vertex plots
1618           meTimeSignalRes_->Fill(recopv.at(ir).recVtx->t() - simpv.at(is).t * simUnit_);
1619           meTimeSignalPull_->Fill((recopv.at(ir).recVtx->t() - simpv.at(is).t * simUnit_) /
1620                                   recopv.at(ir).recVtx->tError());
1621           if (optionalPlots_) {
1622             meRecoPosInSimCollection_->Fill(recopv.at(ir).sim);
1623             meRecoPosInRecoOrigCollection_->Fill(recopv.at(ir).OriginalIndex);
1624           }
1625         }
1626         if (simpv.at(is).eventId.bunchCrossing() == 0 && simpv.at(is).eventId.event() == 0) {
1627           if (!recopv.at(ir).is_signal()) {
1628             edm::LogWarning("Primary4DVertexValidation")
1629                 << "Reco vtx leading match inconsistent: BX/ID " << simpv.at(is).eventId.bunchCrossing() << " "
1630                 << simpv.at(is).eventId.event();
1631           }
1632           meRecoPVPosSignal_->Fill(
1633               recopv.at(ir).OriginalIndex);  // position in reco vtx correction associated to sim signal
1634           if (!signal_is_highest_pt) {
1635             meRecoPVPosSignalNotHighestPt_->Fill(
1636                 recopv.at(ir).OriginalIndex);  // position in reco vtx correction associated to sim signal
1637           }
1638         }
1639 
1640         if (debug_) {
1641           edm::LogPrint("Primary4DVertexValidation") << "*** Matching RECO: " << ir << "with SIM: " << is << " ***";
1642           edm::LogPrint("Primary4DVertexValidation") << "Match Quality is " << recopv.at(ir).matchQuality;
1643           edm::LogPrint("Primary4DVertexValidation") << "****";
1644         }
1645       }
1646     }
1647   }
1648 
1649   //dz histos
1650   for (unsigned int iv = 0; iv < recVtxs->size() - 1; iv++) {
1651     if (recVtxs->at(iv).ndof() > selNdof_) {
1652       double mindistance_realreal = 1e10;
1653 
1654       for (unsigned int jv = iv; jv < recVtxs->size(); jv++) {
1655         if ((!(jv == iv)) && select(recVtxs->at(jv))) {
1656           double dz = recVtxs->at(iv).z() - recVtxs->at(jv).z();
1657           double dtsigma = std::sqrt(recVtxs->at(iv).covariance(3, 3) + recVtxs->at(jv).covariance(3, 3));
1658           double dt = (std::abs(dz) <= deltaZcut_ && dtsigma > 0.)
1659                           ? (recVtxs->at(iv).t() - recVtxs->at(jv).t()) / dtsigma
1660                           : -9999.;
1661           if (recopv.at(iv).is_real() && recopv.at(jv).is_real()) {
1662             meDeltaZrealreal_->Fill(std::abs(dz));
1663             if (dt != -9999.) {
1664               meDeltaTrealreal_->Fill(std::abs(dt));
1665             }
1666             if (std::abs(dz) < std::abs(mindistance_realreal)) {
1667               mindistance_realreal = dz;
1668             }
1669           } else if (recopv.at(iv).is_fake() && recopv.at(jv).is_fake()) {
1670             meDeltaZfakefake_->Fill(std::abs(dz));
1671             if (dt != -9999.) {
1672               meDeltaTfakefake_->Fill(std::abs(dt));
1673             }
1674           }
1675         }
1676       }
1677 
1678       double mindistance_fakereal = 1e10;
1679       double mindistance_realfake = 1e10;
1680       for (unsigned int jv = 0; jv < recVtxs->size(); jv++) {
1681         if ((!(jv == iv)) && select(recVtxs->at(jv))) {
1682           double dz = recVtxs->at(iv).z() - recVtxs->at(jv).z();
1683           double dtsigma = std::sqrt(recVtxs->at(iv).covariance(3, 3) + recVtxs->at(jv).covariance(3, 3));
1684           double dt = (std::abs(dz) <= deltaZcut_ && dtsigma > 0.)
1685                           ? (recVtxs->at(iv).t() - recVtxs->at(jv).t()) / dtsigma
1686                           : -9999.;
1687           if (recopv.at(iv).is_fake() && recopv.at(jv).is_real()) {
1688             meDeltaZfakereal_->Fill(std::abs(dz));
1689             if (dt != -9999.) {
1690               meDeltaTfakereal_->Fill(std::abs(dt));
1691             }
1692             if (std::abs(dz) < std::abs(mindistance_fakereal)) {
1693               mindistance_fakereal = dz;
1694             }
1695           }
1696 
1697           if (recopv.at(iv).is_real() && recopv.at(jv).is_fake() && (std::abs(dz) < std::abs(mindistance_realfake))) {
1698             mindistance_realfake = dz;
1699           }
1700         }
1701       }
1702     }  //ndof
1703   }
1704 
1705 }  // end of analyze
1706 
1707 void Primary4DVertexValidation::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
1708   edm::ParameterSetDescription desc;
1709 
1710   desc.add<std::string>("folder", "MTD/Vertices");
1711   desc.add<edm::InputTag>("TPtoRecoTrackAssoc", edm::InputTag("trackingParticleRecoTrackAsssociation"));
1712   desc.add<edm::InputTag>("mtdTracks", edm::InputTag("trackExtenderWithMTD"));
1713   desc.add<edm::InputTag>("SimTag", edm::InputTag("mix", "MergedTrackTruth"));
1714   desc.add<edm::InputTag>("offlineBS", edm::InputTag("offlineBeamSpot"));
1715   desc.add<edm::InputTag>("offline4DPV", edm::InputTag("offlinePrimaryVertices4D"));
1716   desc.add<edm::InputTag>("trackAssocSrc", edm::InputTag("trackExtenderWithMTD:generalTrackassoc"))
1717       ->setComment("Association between General and MTD Extended tracks");
1718   desc.add<edm::InputTag>("pathLengthSrc", edm::InputTag("trackExtenderWithMTD:generalTrackPathLength"));
1719   desc.add<edm::InputTag>("momentumSrc", edm::InputTag("trackExtenderWithMTD:generalTrackp"));
1720   desc.add<edm::InputTag>("tmtd", edm::InputTag("trackExtenderWithMTD:generalTracktmtd"));
1721   desc.add<edm::InputTag>("timeSrc", edm::InputTag("trackExtenderWithMTD:generalTracktmtd"));
1722   desc.add<edm::InputTag>("sigmaSrc", edm::InputTag("trackExtenderWithMTD:generalTracksigmatmtd"));
1723   desc.add<edm::InputTag>("t0PID", edm::InputTag("tofPID:t0"));
1724   desc.add<edm::InputTag>("t0SafePID", edm::InputTag("tofPID:t0safe"));
1725   desc.add<edm::InputTag>("sigmat0SafePID", edm::InputTag("tofPID:sigmat0safe"));
1726   desc.add<edm::InputTag>("trackMVAQual", edm::InputTag("mtdTrackQualityMVA:mtdQualMVA"));
1727   desc.add<edm::InputTag>("tofPi", edm::InputTag("trackExtenderWithMTD:generalTrackTofPi"));
1728   desc.add<edm::InputTag>("tofK", edm::InputTag("trackExtenderWithMTD:generalTrackTofK"));
1729   desc.add<edm::InputTag>("tofP", edm::InputTag("trackExtenderWithMTD:generalTrackTofP"));
1730   desc.add<edm::InputTag>("probPi", edm::InputTag("tofPID:probPi"));
1731   desc.add<edm::InputTag>("probK", edm::InputTag("tofPID:probK"));
1732   desc.add<edm::InputTag>("probP", edm::InputTag("tofPID:probP"));
1733   desc.add<bool>("useOnlyChargedTracks", true);
1734   desc.addUntracked<bool>("debug", false);
1735   desc.addUntracked<bool>("optionalPlots", false);
1736   desc.add<double>("trackweightTh", 0.5);
1737   desc.add<double>("mvaTh", 0.01);
1738   desc.add<double>("minProbHeavy", 0.75);
1739 
1740   //lineDensity parameters have been obtained by fitting the distribution of the z position of the vertices,
1741   //using a 200k single mu ptGun sample (gaussian fit)
1742   std::vector<double> lDP;
1743   lDP.push_back(1.87);
1744   lDP.push_back(0.);
1745   lDP.push_back(42.5);
1746   desc.add<std::vector<double>>("lineDensityPar", lDP);
1747   descriptions.add("vertices4DValid", desc);
1748 }
1749 
1750 const bool Primary4DVertexValidation::mvaTPSel(const TrackingParticle& tp) {
1751   bool match = false;
1752   if (tp.status() != 1) {
1753     return match;
1754   }
1755   match = tp.charge() != 0 && tp.pt() > pTcut_ && std::abs(tp.eta()) < etacutGEN_;
1756   return match;
1757 }
1758 
1759 const bool Primary4DVertexValidation::mvaRecSel(const reco::TrackBase& trk,
1760                                                 const reco::Vertex& vtx,
1761                                                 const double& t0,
1762                                                 const double& st0) {
1763   bool match = false;
1764   match = trk.pt() > pTcut_ && std::abs(trk.eta()) < etacutREC_ && std::abs(trk.vz() - vtx.z()) <= deltaZcut_;
1765   if (st0 > 0.) {
1766     match = match && std::abs(t0 - vtx.t()) < 3. * st0;
1767   }
1768   return match;
1769 }
1770 
1771 DEFINE_FWK_MODULE(Primary4DVertexValidation);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!