Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:33:35 

0001 #include "Validation/RecoVertex/interface/PrimaryVertexAnalyzer4PUSlimmed.h"
0002 
0003 #include "FWCore/Framework/interface/Event.h"
0004 #include "FWCore/Framework/interface/EventSetup.h"
0005 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007 
0008 // reco track and vertex
0009 #include "DataFormats/VertexReco/interface/Vertex.h"
0010 #include "RecoVertex/VertexPrimitives/interface/TransientVertex.h"
0011 
0012 // TrackingParticle
0013 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h"
0014 #include "SimDataFormats/TrackingAnalysis/interface/TrackingVertexContainer.h"
0015 
0016 // associator
0017 #include "SimTracker/VertexAssociation/interface/calculateVertexSharedTracks.h"
0018 
0019 // DQM
0020 #include "DQMServices/Core/interface/DQMStore.h"
0021 
0022 #include <numeric>
0023 
0024 namespace {
0025   template <typename T, size_t N>
0026   std::array<T, N + 1> makeLogBins(const double min, const double max) {
0027     const double minLog10 = std::log10(min);
0028     const double maxLog10 = std::log10(max);
0029     const double width = (maxLog10 - minLog10) / N;
0030     std::array<T, N + 1> ret;
0031     ret[0] = std::pow(10, minLog10);
0032     const double mult = std::pow(10, width);
0033     for (size_t i = 1; i <= N; ++i) {
0034       ret[i] = ret[i - 1] * mult;
0035     }
0036     return ret;
0037   }
0038 }  // namespace
0039 
0040 //
0041 // constructors and destructor
0042 //
0043 PrimaryVertexAnalyzer4PUSlimmed::PrimaryVertexAnalyzer4PUSlimmed(const edm::ParameterSet& iConfig)
0044     : verbose_(iConfig.getUntrackedParameter<bool>("verbose", false)),
0045       use_only_charged_tracks_(iConfig.getUntrackedParameter<bool>("use_only_charged_tracks", true)),
0046       do_generic_sim_plots_(iConfig.getUntrackedParameter<bool>("do_generic_sim_plots")),
0047       root_folder_(iConfig.getUntrackedParameter<std::string>("root_folder", "Validation/Vertices")),
0048       vecPileupSummaryInfoToken_(consumes<std::vector<PileupSummaryInfo>>(edm::InputTag(std::string("addPileupInfo")))),
0049       trackingParticleCollectionToken_(consumes<TrackingParticleCollection>(
0050           iConfig.getUntrackedParameter<edm::InputTag>("trackingParticleCollection"))),
0051       trackingVertexCollectionToken_(
0052           consumes<TrackingVertexCollection>(iConfig.getUntrackedParameter<edm::InputTag>("trackingVertexCollection"))),
0053       simToRecoAssociationToken_(
0054           consumes<reco::SimToRecoCollection>(iConfig.getUntrackedParameter<edm::InputTag>("trackAssociatorMap"))),
0055       recoToSimAssociationToken_(
0056           consumes<reco::RecoToSimCollection>(iConfig.getUntrackedParameter<edm::InputTag>("trackAssociatorMap"))),
0057       vertexAssociatorToken_(consumes<reco::VertexToTrackingVertexAssociator>(
0058           iConfig.getUntrackedParameter<edm::InputTag>("vertexAssociator"))),
0059       nPUbins_(iConfig.getParameter<unsigned int>("nPUbins")) {
0060   reco_vertex_collections_ = iConfig.getParameter<std::vector<edm::InputTag>>("vertexRecoCollections");
0061   for (auto const& l : reco_vertex_collections_) {
0062     reco_vertex_collection_tokens_.push_back(
0063         edm::EDGetTokenT<edm::View<reco::Vertex>>(consumes<edm::View<reco::Vertex>>(l)));
0064   }
0065   errorPrintedForColl_.resize(reco_vertex_collections_.size(), false);
0066 }
0067 
0068 PrimaryVertexAnalyzer4PUSlimmed::~PrimaryVertexAnalyzer4PUSlimmed() {}
0069 
0070 //
0071 // member functions
0072 //
0073 void PrimaryVertexAnalyzer4PUSlimmed::bookHistograms(DQMStore::IBooker& i,
0074                                                      edm::Run const& iRun,
0075                                                      edm::EventSetup const& iSetup) {
0076   // TODO(rovere) make this booking method shorter and smarter,
0077   // factorizing similar histograms with different prefix in a single
0078   // method call.
0079   float log_bins[31] = {0.0,  0.0002, 0.0004, 0.0006, 0.0008, 0.001, 0.002, 0.004, 0.006, 0.008, 0.01,
0080                         0.02, 0.04,   0.06,   0.08,   0.1,    0.2,   0.4,   0.6,   0.8,   1.0,   2.0,
0081                         4.0,  6.0,    8.0,    10.0,   20.0,   40.0,  60.0,  80.0,  100.0};
0082   auto const log_mergez_bins = makeLogBins<float, 16>(1e-4, 1);  // 4 bins / 10x
0083 
0084   auto const log_pt_bins = makeLogBins<float, 100>(0.1, 1e4);
0085   float log_pt2_bins[16] = {
0086       0.0, 0.1, 0.5, 1.0, 2.0, 5.0, 10.0, 20.0, 50.0, 100.0, 200.0, 500.0, 1000.0, 2000.0, 5000.0, 10000.0};
0087   float log_ntrk_bins[25] = {0.,   2.0,  4.0,  6.0,  8.0,  10.,  12.0, 14.0, 16.0, 18.0,  22.0,  26.0, 30.0,
0088                              35.0, 40.0, 45.0, 50.0, 55.0, 60.0, 70.0, 80.0, 90.0, 100.0, 150.0, 200.0};
0089 
0090   // TODO(rovere) Possibly change or add the main DQMStore booking
0091   // interface to allow booking a TProfile with variable bin-width
0092   // using an array of floats, as done for the TH1F case, not of
0093   // doubles.
0094   double log_pt2_bins_double[16] = {
0095       0.0, 0.1, 0.5, 1.0, 2.0, 5.0, 10.0, 20.0, 50.0, 100.0, 200.0, 500.0, 1000.0, 2000.0, 5000.0, 10000.0};
0096 
0097   i.setCurrentFolder(root_folder_);
0098   if (do_generic_sim_plots_) {
0099     mes_["root_folder"]["GenVtx_vs_BX"] =
0100         i.book2D("GenVtx_vs_BX", "GenVtx_vs_BX", 16, -12.5, 3.5, nPUbins_, 0., double(nPUbins_));
0101     // Generated Primary Vertex Plots
0102     mes_["root_folder"]["GenPV_X"] = i.book1D("GenPV_X", "GeneratedPV_X", 120, -0.6, 0.6);
0103     mes_["root_folder"]["GenPV_Y"] = i.book1D("GenPV_Y", "GeneratedPV_Y", 120, -0.6, 0.6);
0104     mes_["root_folder"]["GenPV_Z"] = i.book1D("GenPV_Z", "GeneratedPV_Z", 120, -60., 60.);
0105     mes_["root_folder"]["GenPV_R"] = i.book1D("GenPV_R", "GeneratedPV_R", 120, 0, 0.6);
0106     mes_["root_folder"]["GenPV_Pt2"] = i.book1D("GenPV_Pt2", "GeneratedPV_Sum-pt2", 15, &log_pt2_bins[0]);
0107     mes_["root_folder"]["GenPV_NumTracks"] =
0108         i.book1D("GenPV_NumTracks", "GeneratedPV_NumTracks", 24, &log_ntrk_bins[0]);
0109     mes_["root_folder"]["GenPV_ClosestDistanceZ"] =
0110         i.book1D("GenPV_ClosestDistanceZ", "GeneratedPV_ClosestDistanceZ", 30, &log_bins[0]);
0111 
0112     // All Generated Vertices, used for efficiency plots
0113     mes_["root_folder"]["GenAllV_NumVertices"] =
0114         i.book1D("GenAllV_NumVertices", "GeneratedAllV_NumVertices", int(nPUbins_ / 2), 0., double(nPUbins_));
0115     mes_["root_folder"]["GenAllV_X"] = i.book1D("GenAllV_X", "GeneratedAllV_X", 120, -0.6, 0.6);
0116     mes_["root_folder"]["GenAllV_Y"] = i.book1D("GenAllV_Y", "GeneratedAllV_Y", 120, -0.6, 0.6);
0117     mes_["root_folder"]["GenAllV_Z"] = i.book1D("GenAllV_Z", "GeneratedAllV_Z", 120, -60, 60);
0118     mes_["root_folder"]["GenAllV_R"] = i.book1D("GenAllV_R", "GeneratedAllV_R", 120, 0, 0.6);
0119     mes_["root_folder"]["GenAllV_Pt2"] = i.book1D("GenAllV_Pt2", "GeneratedAllV_Sum-pt2", 15, &log_pt2_bins[0]);
0120     mes_["root_folder"]["GenAllV_NumTracks"] =
0121         i.book1D("GenAllV_NumTracks", "GeneratedAllV_NumTracks", 24, &log_ntrk_bins[0]);
0122     mes_["root_folder"]["GenAllV_ClosestDistanceZ"] =
0123         i.book1D("GenAllV_ClosestDistanceZ", "GeneratedAllV_ClosestDistanceZ", 30, &log_bins[0]);
0124     mes_["root_folder"]["GenAllV_PairDistanceZ"] =
0125         i.book1D("GenAllV_PairDistanceZ", "GeneratedAllV_PairDistanceZ", 1000, 0, 20);
0126     mes_["root_folder"]["SignalIsHighestPt2"] = i.book1D("SignalIsHighestPt2", "SignalIsHighestPt2", 2, -0.5, 1.5);
0127   }
0128 
0129   for (auto const& l : reco_vertex_collections_) {
0130     std::string label = l.label();
0131     std::string current_folder = root_folder_ + "/" + label;
0132     i.setCurrentFolder(current_folder);
0133 
0134     auto book1d = [&](const char* name, int bins, double min, double max) {
0135       mes_[label][name] = i.book1D(name, name, bins, min, max);
0136     };
0137     auto book1dlogx = [&](const char* name, int bins, const float* xbinedges) {
0138       mes_[label][name] = i.book1D(name, name, bins, xbinedges);
0139     };
0140     auto book2d = [&](const char* name, int xbins, double xmin, double xmax, int ybins, double ymin, double ymax) {
0141       mes_[label][name] = i.book2D(name, name, xbins, xmin, xmax, ybins, ymin, ymax);
0142     };
0143     auto book2dlogx = [&](const char* name, int xbins, const float* xbinedges, int ybins, double ymin, double ymax) {
0144       auto me = i.book2D(name, name, xbins, xbinedges[0], xbinedges[xbins], ybins, ymin, ymax);
0145       me->getTH2F()->GetXaxis()->Set(xbins, xbinedges);
0146       mes_[label][name] = me;
0147     };
0148 
0149     mes_[label]["RecoVtx_vs_GenVtx"] = i.bookProfile(
0150         "RecoVtx_vs_GenVtx", "RecoVtx_vs_GenVtx", nPUbins_, 0., double(nPUbins_), nPUbins_, 0., double(nPUbins_));
0151     mes_[label]["MatchedRecoVtx_vs_GenVtx"] = i.bookProfile("MatchedRecoVtx_vs_GenVtx",
0152                                                             "MatchedRecoVtx_vs_GenVtx",
0153                                                             int(nPUbins_ / 2),
0154                                                             0.,
0155                                                             double(nPUbins_),
0156                                                             nPUbins_,
0157                                                             0.,
0158                                                             double(nPUbins_));
0159     mes_[label]["KindOfSignalPV"] = i.book1D("KindOfSignalPV", "KindOfSignalPV", 9, -0.5, 8.5);
0160     mes_[label]["KindOfSignalPV"]->getTH1()->GetXaxis()->SetBinLabel(1, "!Highest!Assoc2Any");
0161     mes_[label]["KindOfSignalPV"]->getTH1()->GetXaxis()->SetBinLabel(2, "Highest!Assoc2Any");
0162     mes_[label]["KindOfSignalPV"]->getTH1()->GetXaxis()->SetBinLabel(3, "!HighestAssoc2First");
0163     mes_[label]["KindOfSignalPV"]->getTH1()->GetXaxis()->SetBinLabel(4, "HighestAssoc2First");
0164     mes_[label]["KindOfSignalPV"]->getTH1()->GetXaxis()->SetBinLabel(5, "!HighestAssoc2!First");
0165     mes_[label]["KindOfSignalPV"]->getTH1()->GetXaxis()->SetBinLabel(6, "HighestAssoc2!First");
0166     mes_[label]["KindOfSignalPV"]->getTH1()->GetXaxis()->SetBinLabel(7, "!HighestAssoc2First");
0167     mes_[label]["KindOfSignalPV"]->getTH1()->GetXaxis()->SetBinLabel(8, "HighestAssoc2First");
0168     mes_[label]["MisTagRate"] = i.book1D("MisTagRate", "MisTagRate", 2, -0.5, 1.5);
0169     mes_[label]["MisTagRate_vs_PU"] =
0170         i.bookProfile("MisTagRate_vs_PU", "MisTagRate_vs_PU", int(nPUbins_ / 2), 0., double(nPUbins_), 2, 0., 1.);
0171     mes_[label]["MisTagRate_vs_sum-pt2"] =
0172         i.bookProfile("MisTagRate_vs_sum-pt2", "MisTagRate_vs_sum-pt2", 15, &log_pt2_bins_double[0], 2, 0., 1.);
0173     mes_[label]["MisTagRate_vs_Z"] = i.bookProfile("MisTagRate_vs_Z", "MisTagRate_vs_Z", 120, -60., 60., 2, 0., 1.);
0174     mes_[label]["MisTagRate_vs_R"] = i.bookProfile("MisTagRate_vs_R", "MisTagRate_vs_R", 120, 0., 0.6, 2, 0., 1.);
0175     mes_[label]["MisTagRate_vs_NumTracks"] = i.bookProfile(
0176         "MisTagRate_vs_NumTracks", "MisTagRate_vs_NumTracks", int(nPUbins_ / 2), 0., double(nPUbins_), 2, 0., 1.);
0177     mes_[label]["MisTagRateSignalIsHighest"] =
0178         i.book1D("MisTagRateSignalIsHighest", "MisTagRateSignalIsHighest", 2, -0.5, 1.5);
0179     mes_[label]["MisTagRateSignalIsHighest_vs_PU"] = i.bookProfile(
0180         "MisTagRateSignalIsHighest_vs_PU", "MisTagRateSignalIsHighest_vs_PU", nPUbins_, 0., double(nPUbins_), 2, 0., 1.);
0181     mes_[label]["MisTagRateSignalIsHighest_vs_sum-pt2"] = i.bookProfile("MisTagRateSignalIsHighest_vs_sum-pt2",
0182                                                                         "MisTagRateSignalIsHighest_vs_sum-pt2",
0183                                                                         15,
0184                                                                         &log_pt2_bins_double[0],
0185                                                                         2,
0186                                                                         0.,
0187                                                                         1.);
0188     mes_[label]["MisTagRateSignalIsHighest_vs_Z"] =
0189         i.bookProfile("MisTagRateSignalIsHighest_vs_Z", "MisTagRateSignalIsHighest_vs_Z", 120, -60., 60., 2, 0., 1.);
0190     mes_[label]["MisTagRateSignalIsHighest_vs_R"] =
0191         i.bookProfile("MisTagRateSignalIsHighest_vs_R", "MisTagRateSignalIsHighest_vs_R", 120, 0., 0.6, 2, 0., 1.);
0192     mes_[label]["MisTagRateSignalIsHighest_vs_NumTracks"] = i.bookProfile("MisTagRateSignalIsHighest_vs_NumTracks",
0193                                                                           "MisTagRateSignalIsHighest_vs_NumTracks",
0194                                                                           int(nPUbins_ / 2),
0195                                                                           0.,
0196                                                                           double(nPUbins_),
0197                                                                           2,
0198                                                                           0.,
0199                                                                           1.);
0200     mes_[label]["MisTagRateSignalIsNotHighest"] =
0201         i.book1D("MisTagRateSignalIsNotHighest", "MisTagRateSignalIsNotHighest", 2, -0.5, 1.5);
0202     mes_[label]["MisTagRateSignalIsNotHighest_vs_PU"] = i.bookProfile("MisTagRateSignalIsNotHighest_vs_PU",
0203                                                                       "MisTagRateSignalIsNotHighest_vs_PU",
0204                                                                       int(nPUbins_ / 2),
0205                                                                       0.,
0206                                                                       double(nPUbins_),
0207                                                                       2,
0208                                                                       0.,
0209                                                                       1.);
0210     mes_[label]["MisTagRateSignalIsNotHighest_vs_sum-pt2"] = i.bookProfile("MisTagRateSignalIsNotHighest_vs_sum-pt2",
0211                                                                            "MisTagRateSignalIsNotHighest_vs_sum-pt2",
0212                                                                            15,
0213                                                                            &log_pt2_bins_double[0],
0214                                                                            2,
0215                                                                            0.,
0216                                                                            1.);
0217     mes_[label]["MisTagRateSignalIsNotHighest_vs_Z"] = i.bookProfile(
0218         "MisTagRateSignalIsNotHighest_vs_Z", "MisTagRateSignalIsNotHighest_vs_Z", 120, -60., 60., 2, 0., 1.);
0219     mes_[label]["MisTagRateSignalIsNotHighest_vs_R"] = i.bookProfile(
0220         "MisTagRateSignalIsNotHighest_vs_R", "MisTagRateSignalIsNotHighest_vs_R", 120, 0., 0.6, 2, 0., 1.);
0221     mes_[label]["MisTagRateSignalIsNotHighest_vs_NumTracks"] =
0222         i.bookProfile("MisTagRateSignalIsNotHighest_vs_NumTracks",
0223                       "MisTagRateSignalIsNotHighest_vs_NumTracks",
0224                       int(nPUbins_ / 2),
0225                       0.,
0226                       double(nPUbins_),
0227                       2,
0228                       0.,
0229                       1.);
0230     mes_[label]["TruePVLocationIndex"] =
0231         i.book1D("TruePVLocationIndex", "TruePVLocationIndexInRecoVertexCollection", 12, -1.5, 10.5);
0232     mes_[label]["TruePVLocationIndexCumulative"] =
0233         i.book1D("TruePVLocationIndexCumulative", "TruePVLocationIndexInRecoVertexCollectionCumulative", 3, -1.5, 1.5);
0234     mes_[label]["TruePVLocationIndexSignalIsHighest"] =
0235         i.book1D("TruePVLocationIndexSignalIsHighest",
0236                  "TruePVLocationIndexSignalIsHighestInRecoVertexCollection",
0237                  12,
0238                  -1.5,
0239                  10.5);
0240     mes_[label]["TruePVLocationIndexSignalIsNotHighest"] =
0241         i.book1D("TruePVLocationIndexSignalIsNotHighest",
0242                  "TruePVLocationIndexSignalIsNotHighestInRecoVertexCollection",
0243                  12,
0244                  -1.5,
0245                  10.5);
0246     // All Generated Vertices. Used for Efficiency plots We kind of
0247     // duplicate plots here in case we want to perform more detailed
0248     // studies on a selection of generated vertices, not on all of them.
0249     mes_[label]["GenAllAssoc2Reco_NumVertices"] = i.book1D(
0250         "GenAllAssoc2Reco_NumVertices", "GeneratedAllAssoc2Reco_NumVertices", int(nPUbins_ / 2), 0., double(nPUbins_));
0251     mes_[label]["GenAllAssoc2Reco_X"] = i.book1D("GenAllAssoc2Reco_X", "GeneratedAllAssoc2Reco_X", 120, -0.6, 0.6);
0252     mes_[label]["GenAllAssoc2Reco_Y"] = i.book1D("GenAllAssoc2Reco_Y", "GeneratedAllAssoc2Reco_Y", 120, -0.6, 0.6);
0253     mes_[label]["GenAllAssoc2Reco_Z"] = i.book1D("GenAllAssoc2Reco_Z", "GeneratedAllAssoc2Reco_Z", 120, -60, 60);
0254     mes_[label]["GenAllAssoc2Reco_R"] = i.book1D("GenAllAssoc2Reco_R", "GeneratedAllAssoc2Reco_R", 120, 0, 0.6);
0255     mes_[label]["GenAllAssoc2Reco_Pt2"] =
0256         i.book1D("GenAllAssoc2Reco_Pt2", "GeneratedAllAssoc2Reco_Sum-pt2", 15, &log_pt2_bins[0]);
0257     mes_[label]["GenAllAssoc2Reco_NumTracks"] =
0258         i.book1D("GenAllAssoc2Reco_NumTracks", "GeneratedAllAssoc2Reco_NumTracks", 24, &log_ntrk_bins[0]);
0259     mes_[label]["GenAllAssoc2Reco_ClosestDistanceZ"] =
0260         i.book1D("GenAllAssoc2Reco_ClosestDistanceZ", "GeneratedAllAssoc2Reco_ClosestDistanceZ", 30, &log_bins[0]);
0261     book1d("GenPVAssoc2RecoPV_X", 120, -0.6, 0.6);
0262     book1d("GenPVAssoc2RecoPV_Y", 120, -0.6, 0.6);
0263     book1d("GenPVAssoc2RecoPV_Z", 120, -60, 60);
0264 
0265     // All Generated Vertices Matched to a Reconstructed vertex. Used
0266     // for Efficiency plots
0267     mes_[label]["GenAllAssoc2RecoMatched_NumVertices"] = i.book1D("GenAllAssoc2RecoMatched_NumVertices",
0268                                                                   "GeneratedAllAssoc2RecoMatched_NumVertices",
0269                                                                   int(nPUbins_ / 2),
0270                                                                   0.,
0271                                                                   double(nPUbins_));
0272     mes_[label]["GenAllAssoc2RecoMatched_X"] =
0273         i.book1D("GenAllAssoc2RecoMatched_X", "GeneratedAllAssoc2RecoMatched_X", 120, -0.6, 0.6);
0274     mes_[label]["GenAllAssoc2RecoMatched_Y"] =
0275         i.book1D("GenAllAssoc2RecoMatched_Y", "GeneratedAllAssoc2RecoMatched_Y", 120, -0.6, 0.6);
0276     mes_[label]["GenAllAssoc2RecoMatched_Z"] =
0277         i.book1D("GenAllAssoc2RecoMatched_Z", "GeneratedAllAssoc2RecoMatched_Z", 120, -60, 60);
0278     mes_[label]["GenAllAssoc2RecoMatched_R"] =
0279         i.book1D("GenAllAssoc2RecoMatched_R", "GeneratedAllAssoc2RecoMatched_R", 120, 0, 0.6);
0280     mes_[label]["GenAllAssoc2RecoMatched_Pt2"] =
0281         i.book1D("GenAllAssoc2RecoMatched_Pt2", "GeneratedAllAssoc2RecoMatched_Sum-pt2", 15, &log_pt2_bins[0]);
0282     mes_[label]["GenAllAssoc2RecoMatched_NumTracks"] =
0283         i.book1D("GenAllAssoc2RecoMatched_NumTracks", "GeneratedAllAssoc2RecoMatched_NumTracks", 24, &log_ntrk_bins[0]);
0284     mes_[label]["GenAllAssoc2RecoMatched_ClosestDistanceZ"] = i.book1D(
0285         "GenAllAssoc2RecoMatched_ClosestDistanceZ", "GeneratedAllAssoc2RecoMatched_ClosestDistanceZ", 30, &log_bins[0]);
0286     book1d("GenPVAssoc2RecoPVMatched_X", 120, -0.6, 0.6);
0287     book1d("GenPVAssoc2RecoPVMatched_Y", 120, -0.6, 0.6);
0288     book1d("GenPVAssoc2RecoPVMatched_Z", 120, -60, 60);
0289 
0290     // All Generated Vertices Multi-Matched to a Reconstructed vertex. Used
0291     // for Duplicate rate plots
0292     mes_[label]["GenAllAssoc2RecoMultiMatched_NumVertices"] = i.book1D("GenAllAssoc2RecoMultiMatched_NumVertices",
0293                                                                        "GeneratedAllAssoc2RecoMultiMatched_NumVertices",
0294                                                                        int(nPUbins_ / 2),
0295                                                                        0.,
0296                                                                        double(nPUbins_));
0297     mes_[label]["GenAllAssoc2RecoMultiMatched_X"] =
0298         i.book1D("GenAllAssoc2RecoMultiMatched_X", "GeneratedAllAssoc2RecoMultiMatched_X", 120, -0.6, 0.6);
0299     mes_[label]["GenAllAssoc2RecoMultiMatched_Y"] =
0300         i.book1D("GenAllAssoc2RecoMultiMatched_Y", "GeneratedAllAssoc2RecoMultiMatched_Y", 120, -0.6, 0.6);
0301     mes_[label]["GenAllAssoc2RecoMultiMatched_Z"] =
0302         i.book1D("GenAllAssoc2RecoMultiMatched_Z", "GeneratedAllAssoc2RecoMultiMatched_Z", 120, -60, 60);
0303     mes_[label]["GenAllAssoc2RecoMultiMatched_R"] =
0304         i.book1D("GenAllAssoc2RecoMultiMatched_R", "GeneratedAllAssoc2RecoMultiMatched_R", 120, 0, 0.6);
0305     mes_[label]["GenAllAssoc2RecoMultiMatched_Pt2"] = i.book1D(
0306         "GenAllAssoc2RecoMultiMatched_Pt2", "GeneratedAllAssoc2RecoMultiMatched_Sum-pt2", 15, &log_pt2_bins[0]);
0307     mes_[label]["GenAllAssoc2RecoMultiMatched_NumTracks"] = i.book1D("GenAllAssoc2RecoMultiMatched_NumTracks",
0308                                                                      "GeneratedAllAssoc2RecoMultiMatched_NumTracks",
0309                                                                      24,
0310                                                                      &log_ntrk_bins[0]);
0311     mes_[label]["GenAllAssoc2RecoMultiMatched_ClosestDistanceZ"] =
0312         i.book1D("GenAllAssoc2RecoMultiMatched_ClosestDistanceZ",
0313                  "GeneratedAllAssoc2RecoMultiMatched_ClosestDistanceZ",
0314                  30,
0315                  &log_bins[0]);
0316 
0317     // All Reco Vertices. Used for {Fake,Duplicate}-Rate plots
0318     mes_[label]["RecoAllAssoc2Gen_NumVertices"] = i.book1D("RecoAllAssoc2Gen_NumVertices",
0319                                                            "ReconstructedAllAssoc2Gen_NumVertices",
0320                                                            int(nPUbins_ / 2),
0321                                                            0.,
0322                                                            double(nPUbins_));
0323     mes_[label]["RecoAllAssoc2Gen_X"] = i.book1D("RecoAllAssoc2Gen_X", "ReconstructedAllAssoc2Gen_X", 120, -0.6, 0.6);
0324     mes_[label]["RecoAllAssoc2Gen_Y"] = i.book1D("RecoAllAssoc2Gen_Y", "ReconstructedAllAssoc2Gen_Y", 120, -0.6, 0.6);
0325     mes_[label]["RecoAllAssoc2Gen_Z"] = i.book1D("RecoAllAssoc2Gen_Z", "ReconstructedAllAssoc2Gen_Z", 120, -60, 60);
0326     mes_[label]["RecoAllAssoc2Gen_R"] = i.book1D("RecoAllAssoc2Gen_R", "ReconstructedAllAssoc2Gen_R", 120, 0, 0.6);
0327     mes_[label]["RecoAllAssoc2Gen_Pt2"] =
0328         i.book1D("RecoAllAssoc2Gen_Pt2", "ReconstructedAllAssoc2Gen_Sum-pt2", 15, &log_pt2_bins[0]);
0329     mes_[label]["RecoAllAssoc2Gen_Ndof"] =
0330         i.book1D("RecoAllAssoc2Gen_Ndof", "ReconstructedAllAssoc2Gen_Ndof", 120, 0., 240.);
0331     mes_[label]["RecoAllAssoc2Gen_NumTracks"] =
0332         i.book1D("RecoAllAssoc2Gen_NumTracks", "ReconstructedAllAssoc2Gen_NumTracks", 24, &log_ntrk_bins[0]);
0333     mes_[label]["RecoAllAssoc2Gen_PU"] =
0334         i.book1D("RecoAllAssoc2Gen_PU", "ReconstructedAllAssoc2Gen_PU", int(nPUbins_ / 2), 0., double(nPUbins_));
0335     mes_[label]["RecoAllAssoc2Gen_ClosestDistanceZ"] =
0336         i.book1D("RecoAllAssoc2Gen_ClosestDistanceZ", "ReconstructedAllAssoc2Gen_ClosestDistanceZ", 30, &log_bins[0]);
0337     mes_[label]["RecoAllAssoc2GenProperties"] =
0338         i.book1D("RecoAllAssoc2GenProperties", "ReconstructedAllAssoc2Gen_Properties", 8, -0.5, 7.5);
0339     mes_[label]["RecoAllAssoc2Gen_PairDistanceZ"] =
0340         i.book1D("RecoAllAssoc2Gen_PairDistanceZ", "RecoAllAssoc2Gen_PairDistanceZ", 1000, 0, 20);
0341 
0342     // All Reconstructed Vertices Matched to a Generated vertex. Used
0343     // for Fake-Rate plots
0344     mes_[label]["RecoAllAssoc2GenMatched_NumVertices"] = i.book1D("RecoAllAssoc2GenMatched_NumVertices",
0345                                                                   "ReconstructedAllAssoc2GenMatched_NumVertices",
0346                                                                   int(nPUbins_ / 2),
0347                                                                   0.,
0348                                                                   double(nPUbins_));
0349     mes_[label]["RecoAllAssoc2GenMatched_X"] =
0350         i.book1D("RecoAllAssoc2GenMatched_X", "ReconstructedAllAssoc2GenMatched_X", 120, -0.6, 0.6);
0351     mes_[label]["RecoAllAssoc2GenMatched_Y"] =
0352         i.book1D("RecoAllAssoc2GenMatched_Y", "ReconstructedAllAssoc2GenMatched_Y", 120, -0.6, 0.6);
0353     mes_[label]["RecoAllAssoc2GenMatched_Z"] =
0354         i.book1D("RecoAllAssoc2GenMatched_Z", "ReconstructedAllAssoc2GenMatched_Z", 120, -60, 60);
0355     mes_[label]["RecoAllAssoc2GenMatched_R"] =
0356         i.book1D("RecoAllAssoc2GenMatched_R", "ReconstructedAllAssoc2GenMatched_R", 120, 0, 0.6);
0357     mes_[label]["RecoAllAssoc2GenMatched_Pt2"] =
0358         i.book1D("RecoAllAssoc2GenMatched_Pt2", "ReconstructedAllAssoc2GenMatched_Sum-pt2", 15, &log_pt2_bins[0]);
0359     mes_[label]["RecoAllAssoc2GenMatched_Ndof"] =
0360         i.book1D("RecoAllAssoc2GenMatched_Ndof", "ReconstructedAllAssoc2GenMatched_Ndof", 120, 0., 240.);
0361     mes_[label]["RecoAllAssoc2GenMatched_NumTracks"] = i.book1D(
0362         "RecoAllAssoc2GenMatched_NumTracks", "ReconstructedAllAssoc2GenMatched_NumTracks", 24, &log_ntrk_bins[0]);
0363     mes_[label]["RecoAllAssoc2GenMatched_PU"] = i.book1D(
0364         "RecoAllAssoc2GenMatched_PU", "ReconstructedAllAssoc2GenMatched_PU", int(nPUbins_ / 2), 0., double(nPUbins_));
0365     mes_[label]["RecoAllAssoc2GenMatched_ClosestDistanceZ"] =
0366         i.book1D("RecoAllAssoc2GenMatched_ClosestDistanceZ",
0367                  "ReconstructedAllAssoc2GenMatched_ClosestDistanceZ",
0368                  30,
0369                  &log_bins[0]);
0370 
0371     // All Reconstructed Vertices  Multi-Matched to a Generated vertex. Used
0372     // for Merge-Rate plots
0373     mes_[label]["RecoAllAssoc2GenMultiMatched_NumVertices"] =
0374         i.book1D("RecoAllAssoc2GenMultiMatched_NumVertices",
0375                  "ReconstructedAllAssoc2GenMultiMatched_NumVertices",
0376                  int(nPUbins_ / 2),
0377                  0.,
0378                  double(nPUbins_));
0379     mes_[label]["RecoAllAssoc2GenMultiMatched_X"] =
0380         i.book1D("RecoAllAssoc2GenMultiMatched_X", "ReconstructedAllAssoc2GenMultiMatched_X", 120, -0.6, 0.6);
0381     mes_[label]["RecoAllAssoc2GenMultiMatched_Y"] =
0382         i.book1D("RecoAllAssoc2GenMultiMatched_Y", "ReconstructedAllAssoc2GenMultiMatched_Y", 120, -0.6, 0.6);
0383     mes_[label]["RecoAllAssoc2GenMultiMatched_Z"] =
0384         i.book1D("RecoAllAssoc2GenMultiMatched_Z", "ReconstructedAllAssoc2GenMultiMatched_Z", 120, -60, 60);
0385     mes_[label]["RecoAllAssoc2GenMultiMatched_R"] =
0386         i.book1D("RecoAllAssoc2GenMultiMatched_R", "ReconstructedAllAssoc2GenMultiMatched_R", 120, 0, 0.6);
0387     mes_[label]["RecoAllAssoc2GenMultiMatched_Pt2"] = i.book1D(
0388         "RecoAllAssoc2GenMultiMatched_Pt2", "ReconstructedAllAssoc2GenMultiMatched_Sum-pt2", 15, &log_pt2_bins[0]);
0389     mes_[label]["RecoAllAssoc2GenMultiMatched_NumTracks"] = i.book1D("RecoAllAssoc2GenMultiMatched_NumTracks",
0390                                                                      "ReconstructedAllAssoc2GenMultiMatched_NumTracks",
0391                                                                      24,
0392                                                                      &log_ntrk_bins[0]);
0393     mes_[label]["RecoAllAssoc2GenMultiMatched_PU"] = i.book1D("RecoAllAssoc2GenMultiMatched_PU",
0394                                                               "ReconstructedAllAssoc2GenMultiMatched_PU",
0395                                                               int(nPUbins_ / 2),
0396                                                               0.,
0397                                                               double(nPUbins_));
0398     mes_[label]["RecoAllAssoc2GenMultiMatched_ClosestDistanceZ"] =
0399         i.book1D("RecoAllAssoc2GenMultiMatched_ClosestDistanceZ",
0400                  "ReconstructedAllAssoc2GenMultiMatched_ClosestDistanceZ",
0401                  log_mergez_bins.size() - 1,
0402                  &log_mergez_bins[0]);
0403 
0404     // All Reconstructed Vertices Matched to a Multi-Matched Gen
0405     // Vertex. Used for Duplicate rate plots done w.r.t. Reco
0406     // Quantities. We basically want to ask how many times a RecoVTX
0407     // has been reconstructed and associated to a SimulatedVTX that
0408     // has been linked to at least another RecoVTX. In this sense this
0409     // RecoVTX is a duplicate of the same, real GenVTX.
0410     mes_[label]["RecoAllAssoc2MultiMatchedGen_NumVertices"] = i.book1D("RecoAllAssoc2MultiMatchedGen_NumVertices",
0411                                                                        "RecoAllAssoc2MultiMatchedGen_NumVertices",
0412                                                                        int(nPUbins_ / 2),
0413                                                                        0.,
0414                                                                        double(nPUbins_));
0415     mes_[label]["RecoAllAssoc2MultiMatchedGen_X"] =
0416         i.book1D("RecoAllAssoc2MultiMatchedGen_X", "RecoAllAssoc2MultiMatchedGen_X", 120, -0.6, 0.6);
0417     mes_[label]["RecoAllAssoc2MultiMatchedGen_Y"] =
0418         i.book1D("RecoAllAssoc2MultiMatchedGen_Y", "RecoAllAssoc2MultiMatchedGen_Y", 120, -0.6, 0.6);
0419     mes_[label]["RecoAllAssoc2MultiMatchedGen_Z"] =
0420         i.book1D("RecoAllAssoc2MultiMatchedGen_Z", "RecoAllAssoc2MultiMatchedGen_Z", 120, -60, 60);
0421     mes_[label]["RecoAllAssoc2MultiMatchedGen_R"] =
0422         i.book1D("RecoAllAssoc2MultiMatchedGen_R", "RecoAllAssoc2MultiMatchedGen_R", 120, 0, 0.6);
0423     mes_[label]["RecoAllAssoc2MultiMatchedGen_Pt2"] =
0424         i.book1D("RecoAllAssoc2MultiMatchedGen_Pt2", "RecoAllAssoc2MultiMatchedGen_Sum-pt2", 15, &log_pt2_bins[0]);
0425     mes_[label]["RecoAllAssoc2MultiMatchedGen_NumTracks"] = i.book1D(
0426         "RecoAllAssoc2MultiMatchedGen_NumTracks", "RecoAllAssoc2MultiMatchedGen_NumTracks", 24, &log_ntrk_bins[0]);
0427     mes_[label]["RecoAllAssoc2MultiMatchedGen_PU"] = i.book1D(
0428         "RecoAllAssoc2MultiMatchedGen_PU", "RecoAllAssoc2MultiMatchedGen_PU", int(nPUbins_ / 2), 0., double(nPUbins_));
0429     mes_[label]["RecoAllAssoc2MultiMatchedGen_ClosestDistanceZ"] =
0430         i.book1D("RecoAllAssoc2MultiMatchedGen_ClosestDistanceZ",
0431                  "RecoAllAssoc2MultiMatchedGen_ClosestDistanceZ",
0432                  30,
0433                  &log_bins[0]);
0434     mes_[label]["RecoAllAssoc2GenSimForMerge_ClosestDistanceZ"] =
0435         i.book1D("RecoAllAssoc2GenSimForMerge_ClosestDistanceZ",
0436                  "RecoAllAssoc2GenSimForMerge_ClosestDistanceZ",
0437                  log_mergez_bins.size() - 1,
0438                  &log_mergez_bins[0]);
0439 
0440     // Resolution and pull histograms
0441     const double resolpt2 = 10;
0442 
0443     const double minPull = -10;
0444     const double maxPull = 10;
0445     const double nPull = 100;
0446 
0447     auto bookResolPull = [&](const std::string& prefix, const double resolx, const double resoly, const double resolz) {
0448       book1d((prefix + "_ResolX").c_str(), 100, -resolx, resolx);
0449       book1d((prefix + "_ResolY").c_str(), 100, -resoly, resoly);
0450       book1d((prefix + "_ResolZ").c_str(), 100, -resolz, resolz);
0451       book1d((prefix + "_ResolPt2").c_str(), 100, -resolpt2, resolpt2);
0452       book1d((prefix + "_PullX").c_str(), 250, -25, 25);
0453       book1d((prefix + "_PullY").c_str(), 250, -25, 25);
0454       book1d((prefix + "_PullZ").c_str(), 250, -25, 25);
0455 
0456       book2d((prefix + "_ResolX_vs_PU").c_str(), int(nPUbins_ / 2), 0., nPUbins_, 100, -resolx, resolx);
0457       book2d((prefix + "_ResolY_vs_PU").c_str(), int(nPUbins_ / 2), 0., double(nPUbins_), 100, -resoly, resoly);
0458       book2d((prefix + "_ResolZ_vs_PU").c_str(), int(nPUbins_ / 2), 0., double(nPUbins_), 100, -resolz, resolz);
0459       book2d((prefix + "_ResolPt2_vs_PU").c_str(), int(nPUbins_ / 2), 0., double(nPUbins_), 100, -resolpt2, resolpt2);
0460       book2d((prefix + "_PullX_vs_PU").c_str(), int(nPUbins_ / 2), 0., double(nPUbins_), nPull, minPull, maxPull);
0461       book2d((prefix + "_PullY_vs_PU").c_str(), int(nPUbins_ / 2), 0., double(nPUbins_), nPull, minPull, maxPull);
0462       book2d((prefix + "_PullZ_vs_PU").c_str(), int(nPUbins_ / 2), 0., double(nPUbins_), nPull, minPull, maxPull);
0463 
0464       book2dlogx((prefix + "_ResolX_vs_NumTracks").c_str(), 24, &log_ntrk_bins[0], 100, -resolx, resolx);
0465       book2dlogx((prefix + "_ResolY_vs_NumTracks").c_str(), 24, &log_ntrk_bins[0], 100, -resoly, resoly);
0466       book2dlogx((prefix + "_ResolZ_vs_NumTracks").c_str(), 24, &log_ntrk_bins[0], 100, -resolz, resolz);
0467       book2dlogx((prefix + "_ResolPt2_vs_NumTracks").c_str(), 24, &log_ntrk_bins[0], 100, -resolpt2, resolpt2);
0468       book2dlogx((prefix + "_PullX_vs_NumTracks").c_str(), 24, &log_ntrk_bins[0], nPull, minPull, maxPull);
0469       book2dlogx((prefix + "_PullY_vs_NumTracks").c_str(), 24, &log_ntrk_bins[0], nPull, minPull, maxPull);
0470       book2dlogx((prefix + "_PullZ_vs_NumTracks").c_str(), 24, &log_ntrk_bins[0], nPull, minPull, maxPull);
0471 
0472       book2d((prefix + "_ResolX_vs_Z").c_str(), 120, -60, 60, 100, -resolx, resolx);
0473       book2d((prefix + "_ResolY_vs_Z").c_str(), 120, -60, 60, 100, -resoly, resoly);
0474       book2d((prefix + "_ResolZ_vs_Z").c_str(), 120, -60, 60, 100, -resolz, resolz);
0475       book2d((prefix + "_ResolPt2_vs_Z").c_str(), 120, -60, 60, 100, -resolpt2, resolpt2);
0476       book2d((prefix + "_PullX_vs_Z").c_str(), 120, -60, 60, nPull, minPull, maxPull);
0477       book2d((prefix + "_PullY_vs_Z").c_str(), 120, -60, 60, nPull, minPull, maxPull);
0478       book2d((prefix + "_PullZ_vs_Z").c_str(), 120, -60, 60, nPull, minPull, maxPull);
0479 
0480       book2dlogx((prefix + "_ResolX_vs_Pt").c_str(), log_pt_bins.size() - 1, &log_pt_bins[0], 100, -resolx, resolx);
0481       book2dlogx((prefix + "_ResolY_vs_Pt").c_str(), log_pt_bins.size() - 1, &log_pt_bins[0], 100, -resoly, resoly);
0482       book2dlogx((prefix + "_ResolZ_vs_Pt").c_str(), log_pt_bins.size() - 1, &log_pt_bins[0], 100, -resolz, resolz);
0483       book2dlogx(
0484           (prefix + "_ResolPt2_vs_Pt").c_str(), log_pt_bins.size() - 1, &log_pt_bins[0], 100, -resolpt2, resolpt2);
0485       book2dlogx((prefix + "_PullX_vs_Pt").c_str(), log_pt_bins.size() - 1, &log_pt_bins[0], nPull, minPull, maxPull);
0486       book2dlogx((prefix + "_PullY_vs_Pt").c_str(), log_pt_bins.size() - 1, &log_pt_bins[0], nPull, minPull, maxPull);
0487       book2dlogx((prefix + "_PullZ_vs_Pt").c_str(), log_pt_bins.size() - 1, &log_pt_bins[0], nPull, minPull, maxPull);
0488     };
0489 
0490     bookResolPull("RecoAllAssoc2GenMatched", 0.1, 0.1, 0.1);        // Non-merged vertices
0491     bookResolPull("RecoAllAssoc2GenMatchedMerged", 0.1, 0.1, 0.1);  // Merged vertices
0492     bookResolPull(
0493         "RecoPVAssoc2GenPVMatched", 0.01, 0.01, 0.01);  // PV, when correctly matched, regardless if merged or not
0494 
0495     // Purity histograms
0496     // Reco PV (vtx0) matched to hard-scatter gen vertex
0497     book1d("RecoPVAssoc2GenPVMatched_Purity", 50, 0, 1);
0498     book1d("RecoPVAssoc2GenPVMatched_Missing", 50, 0, 1);
0499     book2d("RecoPVAssoc2GenPVMatched_Purity_vs_Index", 100, 0, 100, 50, 0, 1);
0500 
0501     // RECO PV (vtx0) not matched to hard-scatter gen vertex
0502     book1d("RecoPVAssoc2GenPVNotMatched_Purity", 50, 0, 1);
0503     book1d("RecoPVAssoc2GenPVNotMatched_Missing", 50, 0, 1);
0504     book2d("RecoPVAssoc2GenPVNotMatched_Purity_vs_Index", 100, 0, 100, 50, 0, 1);
0505 
0506     // Purity vs. fake rate
0507     book1d("RecoAllAssoc2Gen_Purity", 50, 0, 1);         // denominator
0508     book1d("RecoAllAssoc2GenMatched_Purity", 50, 0, 1);  // 1-numerator
0509 
0510     // Vertex sum(pt2)
0511     // The first two are orthogonal (i.e. their sum includes all reco vertices)
0512     book1dlogx("RecoAssoc2GenPVMatched_Pt2", 15, &log_pt2_bins[0]);
0513     book1dlogx("RecoAssoc2GenPVNotMatched_Pt2", 15, &log_pt2_bins[0]);
0514 
0515     book1dlogx("RecoAssoc2GenPVMatchedNotHighest_Pt2", 15, &log_pt2_bins[0]);
0516     book1dlogx("RecoAssoc2GenPVNotMatched_GenPVTracksRemoved_Pt2", 15, &log_pt2_bins[0]);
0517 
0518     // Shared tracks
0519     book1d("RecoAllAssoc2GenSingleMatched_SharedTrackFractionReco", 50, 0, 1);
0520     book1d("RecoAllAssoc2GenMultiMatched_SharedTrackFractionReco", 50, 0, 1);
0521     book1d("RecoAllAssoc2GenSingleMatched_SharedTrackFractionRecoMatched", 50, 0, 1);
0522     book1d("RecoAllAssoc2GenMultiMatched_SharedTrackFractionRecoMatched", 50, 0, 1);
0523 
0524     book1d("RecoAllAssoc2GenSingleMatched_SharedTrackFractionSim", 50, 0, 1);
0525     book1d("RecoAllAssoc2GenMultiMatched_SharedTrackFractionSim", 50, 0, 1);
0526     book1d("RecoAllAssoc2GenSingleMatched_SharedTrackFractionSimMatched", 50, 0, 1);
0527     book1d("RecoAllAssoc2GenMultiMatched_SharedTrackFractionSimMatched", 50, 0, 1);
0528   }
0529 }
0530 
0531 void PrimaryVertexAnalyzer4PUSlimmed::fillGenericGenVertexHistograms(const simPrimaryVertex& v) {
0532   if (v.eventId.event() == 0) {
0533     mes_["root_folder"]["GenPV_X"]->Fill(v.x);
0534     mes_["root_folder"]["GenPV_Y"]->Fill(v.y);
0535     mes_["root_folder"]["GenPV_Z"]->Fill(v.z);
0536     mes_["root_folder"]["GenPV_R"]->Fill(v.r);
0537     mes_["root_folder"]["GenPV_Pt2"]->Fill(v.ptsq);
0538     mes_["root_folder"]["GenPV_NumTracks"]->Fill(v.nGenTrk);
0539     if (v.closest_vertex_distance_z > 0.)
0540       mes_["root_folder"]["GenPV_ClosestDistanceZ"]->Fill(v.closest_vertex_distance_z);
0541   }
0542   mes_["root_folder"]["GenAllV_X"]->Fill(v.x);
0543   mes_["root_folder"]["GenAllV_Y"]->Fill(v.y);
0544   mes_["root_folder"]["GenAllV_Z"]->Fill(v.z);
0545   mes_["root_folder"]["GenAllV_R"]->Fill(v.r);
0546   mes_["root_folder"]["GenAllV_Pt2"]->Fill(v.ptsq);
0547   mes_["root_folder"]["GenAllV_NumTracks"]->Fill(v.nGenTrk);
0548   if (v.closest_vertex_distance_z > 0.)
0549     mes_["root_folder"]["GenAllV_ClosestDistanceZ"]->Fill(v.closest_vertex_distance_z);
0550 }
0551 
0552 void PrimaryVertexAnalyzer4PUSlimmed::fillRecoAssociatedGenVertexHistograms(
0553     const std::string& label, const PrimaryVertexAnalyzer4PUSlimmed::simPrimaryVertex& v) {
0554   mes_[label]["GenAllAssoc2Reco_X"]->Fill(v.x);
0555   mes_[label]["GenAllAssoc2Reco_Y"]->Fill(v.y);
0556   mes_[label]["GenAllAssoc2Reco_Z"]->Fill(v.z);
0557   mes_[label]["GenAllAssoc2Reco_R"]->Fill(v.r);
0558   mes_[label]["GenAllAssoc2Reco_Pt2"]->Fill(v.ptsq);
0559   mes_[label]["GenAllAssoc2Reco_NumTracks"]->Fill(v.nGenTrk);
0560   if (v.closest_vertex_distance_z > 0.)
0561     mes_[label]["GenAllAssoc2Reco_ClosestDistanceZ"]->Fill(v.closest_vertex_distance_z);
0562   if (!v.rec_vertices.empty()) {
0563     mes_[label]["GenAllAssoc2RecoMatched_X"]->Fill(v.x);
0564     mes_[label]["GenAllAssoc2RecoMatched_Y"]->Fill(v.y);
0565     mes_[label]["GenAllAssoc2RecoMatched_Z"]->Fill(v.z);
0566     mes_[label]["GenAllAssoc2RecoMatched_R"]->Fill(v.r);
0567     mes_[label]["GenAllAssoc2RecoMatched_Pt2"]->Fill(v.ptsq);
0568     mes_[label]["GenAllAssoc2RecoMatched_NumTracks"]->Fill(v.nGenTrk);
0569     if (v.closest_vertex_distance_z > 0.)
0570       mes_[label]["GenAllAssoc2RecoMatched_ClosestDistanceZ"]->Fill(v.closest_vertex_distance_z);
0571   }
0572   if (v.rec_vertices.size() > 1) {
0573     mes_[label]["GenAllAssoc2RecoMultiMatched_X"]->Fill(v.x);
0574     mes_[label]["GenAllAssoc2RecoMultiMatched_Y"]->Fill(v.y);
0575     mes_[label]["GenAllAssoc2RecoMultiMatched_Z"]->Fill(v.z);
0576     mes_[label]["GenAllAssoc2RecoMultiMatched_R"]->Fill(v.r);
0577     mes_[label]["GenAllAssoc2RecoMultiMatched_Pt2"]->Fill(v.ptsq);
0578     mes_[label]["GenAllAssoc2RecoMultiMatched_NumTracks"]->Fill(v.nGenTrk);
0579     if (v.closest_vertex_distance_z > 0.)
0580       mes_[label]["GenAllAssoc2RecoMultiMatched_ClosestDistanceZ"]->Fill(v.closest_vertex_distance_z);
0581   }
0582 }
0583 
0584 void PrimaryVertexAnalyzer4PUSlimmed::fillRecoAssociatedGenPVHistograms(
0585     const std::string& label, const PrimaryVertexAnalyzer4PUSlimmed::simPrimaryVertex& v, bool genPVMatchedToRecoPV) {
0586   mes_[label]["GenPVAssoc2RecoPV_X"]->Fill(v.x);
0587   mes_[label]["GenPVAssoc2RecoPV_Y"]->Fill(v.y);
0588   mes_[label]["GenPVAssoc2RecoPV_Z"]->Fill(v.z);
0589   if (genPVMatchedToRecoPV) {
0590     mes_[label]["GenPVAssoc2RecoPVMatched_X"]->Fill(v.x);
0591     mes_[label]["GenPVAssoc2RecoPVMatched_Y"]->Fill(v.y);
0592     mes_[label]["GenPVAssoc2RecoPVMatched_Z"]->Fill(v.z);
0593   }
0594 }
0595 
0596 void PrimaryVertexAnalyzer4PUSlimmed::fillGenAssociatedRecoVertexHistograms(
0597     const std::string& label, int num_pileup_vertices, PrimaryVertexAnalyzer4PUSlimmed::recoPrimaryVertex& v) {
0598   mes_[label]["RecoAllAssoc2Gen_X"]->Fill(v.x);
0599   mes_[label]["RecoAllAssoc2Gen_Y"]->Fill(v.y);
0600   mes_[label]["RecoAllAssoc2Gen_Z"]->Fill(v.z);
0601   mes_[label]["RecoAllAssoc2Gen_R"]->Fill(v.r);
0602   mes_[label]["RecoAllAssoc2Gen_Pt2"]->Fill(v.ptsq);
0603   mes_[label]["RecoAllAssoc2Gen_Ndof"]->Fill(v.recVtx->ndof());
0604   mes_[label]["RecoAllAssoc2Gen_NumTracks"]->Fill(v.nRecoTrk);
0605   mes_[label]["RecoAllAssoc2Gen_PU"]->Fill(num_pileup_vertices);
0606   mes_[label]["RecoAllAssoc2Gen_Purity"]->Fill(v.purity);
0607   if (v.closest_vertex_distance_z > 0.)
0608     mes_[label]["RecoAllAssoc2Gen_ClosestDistanceZ"]->Fill(v.closest_vertex_distance_z);
0609   if (!v.sim_vertices.empty()) {
0610     v.kind_of_vertex |= recoPrimaryVertex::MATCHED;
0611     mes_[label]["RecoAllAssoc2GenMatched_X"]->Fill(v.x);
0612     mes_[label]["RecoAllAssoc2GenMatched_Y"]->Fill(v.y);
0613     mes_[label]["RecoAllAssoc2GenMatched_Z"]->Fill(v.z);
0614     mes_[label]["RecoAllAssoc2GenMatched_R"]->Fill(v.r);
0615     mes_[label]["RecoAllAssoc2GenMatched_Pt2"]->Fill(v.ptsq);
0616     mes_[label]["RecoAllAssoc2GenMatched_Ndof"]->Fill(v.recVtx->ndof());
0617     mes_[label]["RecoAllAssoc2GenMatched_NumTracks"]->Fill(v.nRecoTrk);
0618     mes_[label]["RecoAllAssoc2GenMatched_PU"]->Fill(num_pileup_vertices);
0619     mes_[label]["RecoAllAssoc2GenMatched_Purity"]->Fill(v.purity);
0620     if (v.closest_vertex_distance_z > 0.)
0621       mes_[label]["RecoAllAssoc2GenMatched_ClosestDistanceZ"]->Fill(v.closest_vertex_distance_z);
0622 
0623     // Fill resolution and pull plots here (as in MultiTrackValidator)
0624     fillResolutionAndPullHistograms(label, num_pileup_vertices, v, false);
0625 
0626     // Now keep track of all RecoVTX associated to a SimVTX that
0627     // itself is associated to more than one RecoVTX, for
0628     // duplicate-rate plots on reco quantities.
0629     if (v.sim_vertices_internal[0]->rec_vertices.size() > 1) {
0630       v.kind_of_vertex |= recoPrimaryVertex::DUPLICATE;
0631       mes_[label]["RecoAllAssoc2MultiMatchedGen_X"]->Fill(v.x);
0632       mes_[label]["RecoAllAssoc2MultiMatchedGen_Y"]->Fill(v.y);
0633       mes_[label]["RecoAllAssoc2MultiMatchedGen_Z"]->Fill(v.z);
0634       mes_[label]["RecoAllAssoc2MultiMatchedGen_R"]->Fill(v.r);
0635       mes_[label]["RecoAllAssoc2MultiMatchedGen_Pt2"]->Fill(v.ptsq);
0636       mes_[label]["RecoAllAssoc2MultiMatchedGen_NumTracks"]->Fill(v.nRecoTrk);
0637       mes_[label]["RecoAllAssoc2MultiMatchedGen_PU"]->Fill(num_pileup_vertices);
0638       if (v.closest_vertex_distance_z > 0.)
0639         mes_[label]["RecoAllAssoc2MultiMatchedGen_ClosestDistanceZ"]->Fill(v.closest_vertex_distance_z);
0640     }
0641     // This is meant to be used as "denominator" for the merge-rate
0642     // plots produced starting from reco quantities. We   enter here
0643     // only if the reco vertex has been associated, since we need info
0644     // from the SimVTX associated to it. In this regard, the final
0645     // merge-rate plot coming from reco is not to be intended as a
0646     // pure efficiency-like plot, since the normalization is biased.
0647     if (v.sim_vertices_internal[0]->closest_vertex_distance_z > 0.)
0648       mes_[label]["RecoAllAssoc2GenSimForMerge_ClosestDistanceZ"]->Fill(
0649           v.sim_vertices_internal[0]->closest_vertex_distance_z);
0650   }
0651   // this plots are meant to be used to compute the merge rate
0652   if (v.sim_vertices.size() > 1) {
0653     v.kind_of_vertex |= recoPrimaryVertex::MERGED;
0654     mes_[label]["RecoAllAssoc2GenMultiMatched_X"]->Fill(v.x);
0655     mes_[label]["RecoAllAssoc2GenMultiMatched_Y"]->Fill(v.y);
0656     mes_[label]["RecoAllAssoc2GenMultiMatched_Z"]->Fill(v.z);
0657     mes_[label]["RecoAllAssoc2GenMultiMatched_R"]->Fill(v.r);
0658     mes_[label]["RecoAllAssoc2GenMultiMatched_Pt2"]->Fill(v.ptsq);
0659     mes_[label]["RecoAllAssoc2GenMultiMatched_NumTracks"]->Fill(v.nRecoTrk);
0660     mes_[label]["RecoAllAssoc2GenMultiMatched_PU"]->Fill(num_pileup_vertices);
0661     if (v.sim_vertices_internal[0]->closest_vertex_distance_z > 0.)
0662       mes_[label]["RecoAllAssoc2GenMultiMatched_ClosestDistanceZ"]->Fill(
0663           v.sim_vertices_internal[0]->closest_vertex_distance_z);
0664   }
0665   mes_[label]["RecoAllAssoc2GenProperties"]->Fill(v.kind_of_vertex);
0666 
0667   std::string prefix;
0668   if (v.sim_vertices.size() == 1) {
0669     prefix = "RecoAllAssoc2GenSingleMatched_SharedTrackFraction";
0670   } else if (v.sim_vertices.size() > 1) {
0671     prefix = "RecoAllAssoc2GenMultiMatched_SharedTrackFraction";
0672   }
0673 
0674   for (size_t i = 0; i < v.sim_vertices.size(); ++i) {
0675     const double sharedTracks = v.sim_vertices_num_shared_tracks[i];
0676     const simPrimaryVertex* simV = v.sim_vertices_internal[i];
0677     mes_[label][prefix + "Reco"]->Fill(sharedTracks / v.nRecoTrk);
0678     mes_[label][prefix + "RecoMatched"]->Fill(sharedTracks / v.num_matched_sim_tracks);
0679     mes_[label][prefix + "Sim"]->Fill(sharedTracks / simV->nGenTrk);
0680     mes_[label][prefix + "SimMatched"]->Fill(sharedTracks / simV->num_matched_reco_tracks);
0681   }
0682 }
0683 
0684 void PrimaryVertexAnalyzer4PUSlimmed::fillResolutionAndPullHistograms(
0685     const std::string& label,
0686     int num_pileup_vertices,
0687     PrimaryVertexAnalyzer4PUSlimmed::recoPrimaryVertex& v,
0688     bool isPV) {
0689   std::string prefix = "RecoAllAssoc2GenMatched";
0690   if (v.sim_vertices_internal.size() > 1) {
0691     prefix += "Merged";
0692   }
0693   if (isPV) {
0694     prefix = "RecoPVAssoc2GenPVMatched";
0695   }
0696 
0697   // Use the best match as defined by the vertex truth associator
0698   // reco-tracks as the best match
0699   const simPrimaryVertex& bestMatch = *(v.sim_vertices_internal[0]);
0700   const double xres = v.x - bestMatch.x;
0701   const double yres = v.y - bestMatch.y;
0702   const double zres = v.z - bestMatch.z;
0703   const double pt2res = v.ptsq - bestMatch.ptsq;
0704 
0705   const double xresol = xres;
0706   const double yresol = yres;
0707   const double zresol = zres;
0708   const double pt2resol = pt2res / v.ptsq;
0709   const double xpull = xres / v.recVtx->xError();
0710   const double ypull = yres / v.recVtx->yError();
0711   const double zpull = zres / v.recVtx->zError();
0712 
0713   mes_[label][prefix + "_ResolX"]->Fill(xresol);
0714   mes_[label][prefix + "_ResolY"]->Fill(yresol);
0715   mes_[label][prefix + "_ResolZ"]->Fill(zresol);
0716   mes_[label][prefix + "_ResolPt2"]->Fill(pt2resol);
0717   mes_[label][prefix + "_PullX"]->Fill(xpull);
0718   mes_[label][prefix + "_PullY"]->Fill(ypull);
0719   mes_[label][prefix + "_PullZ"]->Fill(zpull);
0720 
0721   mes_[label][prefix + "_ResolX_vs_PU"]->Fill(num_pileup_vertices, xresol);
0722   mes_[label][prefix + "_ResolY_vs_PU"]->Fill(num_pileup_vertices, yresol);
0723   mes_[label][prefix + "_ResolZ_vs_PU"]->Fill(num_pileup_vertices, zresol);
0724   mes_[label][prefix + "_ResolPt2_vs_PU"]->Fill(num_pileup_vertices, pt2resol);
0725   mes_[label][prefix + "_PullX_vs_PU"]->Fill(num_pileup_vertices, xpull);
0726   mes_[label][prefix + "_PullY_vs_PU"]->Fill(num_pileup_vertices, ypull);
0727   mes_[label][prefix + "_PullZ_vs_PU"]->Fill(num_pileup_vertices, zpull);
0728 
0729   mes_[label][prefix + "_ResolX_vs_NumTracks"]->Fill(v.nRecoTrk, xresol);
0730   mes_[label][prefix + "_ResolY_vs_NumTracks"]->Fill(v.nRecoTrk, yresol);
0731   mes_[label][prefix + "_ResolZ_vs_NumTracks"]->Fill(v.nRecoTrk, zresol);
0732   mes_[label][prefix + "_ResolPt2_vs_NumTracks"]->Fill(v.nRecoTrk, pt2resol);
0733   mes_[label][prefix + "_PullX_vs_NumTracks"]->Fill(v.nRecoTrk, xpull);
0734   mes_[label][prefix + "_PullY_vs_NumTracks"]->Fill(v.nRecoTrk, ypull);
0735   mes_[label][prefix + "_PullZ_vs_NumTracks"]->Fill(v.nRecoTrk, zpull);
0736 
0737   mes_[label][prefix + "_ResolX_vs_Z"]->Fill(v.z, xresol);
0738   mes_[label][prefix + "_ResolY_vs_Z"]->Fill(v.z, yresol);
0739   mes_[label][prefix + "_ResolZ_vs_Z"]->Fill(v.z, zresol);
0740   mes_[label][prefix + "_ResolPt2_vs_Z"]->Fill(v.z, pt2resol);
0741   mes_[label][prefix + "_PullX_vs_Z"]->Fill(v.z, xpull);
0742   mes_[label][prefix + "_PullY_vs_Z"]->Fill(v.z, ypull);
0743   mes_[label][prefix + "_PullZ_vs_Z"]->Fill(v.z, zpull);
0744 
0745   mes_[label][prefix + "_ResolX_vs_Pt"]->Fill(v.pt, xresol);
0746   mes_[label][prefix + "_ResolY_vs_Pt"]->Fill(v.pt, yresol);
0747   mes_[label][prefix + "_ResolZ_vs_Pt"]->Fill(v.pt, zresol);
0748   mes_[label][prefix + "_ResolPt2_vs_Pt"]->Fill(v.pt, pt2resol);
0749   mes_[label][prefix + "_PullX_vs_Pt"]->Fill(v.pt, xpull);
0750   mes_[label][prefix + "_PullY_vs_Pt"]->Fill(v.pt, ypull);
0751   mes_[label][prefix + "_PullZ_vs_Pt"]->Fill(v.pt, zpull);
0752 }
0753 
0754 bool PrimaryVertexAnalyzer4PUSlimmed::matchRecoTrack2SimSignal(const reco::TrackBaseRef& recoTrack) {
0755   auto found = r2s_->find(recoTrack);
0756 
0757   // reco track not matched to any TP
0758   if (found == r2s_->end())
0759     return false;
0760 
0761   // reco track matched to some TP from signal vertex
0762   for (const auto& tp : found->val) {
0763     if (tp.first->eventId().bunchCrossing() == 0 && tp.first->eventId().event() == 0)
0764       return true;
0765   }
0766 
0767   // reco track not matched to any TP from signal vertex
0768   return false;
0769 }
0770 
0771 void PrimaryVertexAnalyzer4PUSlimmed::calculatePurityAndFillHistograms(const std::string& label,
0772                                                                        std::vector<recoPrimaryVertex>& recopvs,
0773                                                                        int genpv_position_in_reco_collection,
0774                                                                        bool signal_is_highest_pt) {
0775   if (recopvs.empty())
0776     return;
0777 
0778   std::vector<double> vtx_sumpt_sigmatched;
0779   std::vector<double> vtx_sumpt2_sigmatched;
0780 
0781   vtx_sumpt_sigmatched.reserve(recopvs.size());
0782   vtx_sumpt2_sigmatched.reserve(recopvs.size());
0783 
0784   // Calculate purity
0785   for (auto& v : recopvs) {
0786     double sumpt_all = 0;
0787     double sumpt_sigmatched = 0;
0788     double sumpt2_sigmatched = 0;
0789     const reco::Vertex* vertex = v.recVtx;
0790     for (auto iTrack = vertex->tracks_begin(); iTrack != vertex->tracks_end(); ++iTrack) {
0791       double pt = (*iTrack)->pt();
0792       sumpt_all += pt;
0793       if (matchRecoTrack2SimSignal(*iTrack)) {
0794         sumpt_sigmatched += pt;
0795         sumpt2_sigmatched += pt * pt;
0796       }
0797     }
0798     v.purity = sumpt_sigmatched / sumpt_all;
0799 
0800     vtx_sumpt_sigmatched.push_back(sumpt_sigmatched);
0801     vtx_sumpt2_sigmatched.push_back(sumpt2_sigmatched);
0802   }
0803 
0804   const double vtxAll_sumpt_sigmatched = std::accumulate(vtx_sumpt_sigmatched.begin(), vtx_sumpt_sigmatched.end(), 0.0);
0805   const double vtxNot0_sumpt_sigmatched = vtxAll_sumpt_sigmatched - vtx_sumpt_sigmatched[0];
0806   const double missing = vtxNot0_sumpt_sigmatched / vtxAll_sumpt_sigmatched;
0807 
0808   // Fill purity
0809   std::string prefix = "RecoPVAssoc2GenPVNotMatched_";
0810   if (genpv_position_in_reco_collection == 0)
0811     prefix = "RecoPVAssoc2GenPVMatched_";
0812 
0813   mes_[label][prefix + "Purity"]->Fill(recopvs[0].purity);
0814   mes_[label][prefix + "Missing"]->Fill(missing);
0815   auto hpurity = mes_[label][prefix + "Purity_vs_Index"];
0816   for (size_t i = 0; i < recopvs.size(); ++i) {
0817     hpurity->Fill(i, recopvs[i].purity);
0818   }
0819 
0820   // Fill sumpt2
0821   for (size_t i = 0; i < recopvs.size(); ++i) {
0822     if (static_cast<int>(i) == genpv_position_in_reco_collection) {
0823       mes_[label]["RecoAssoc2GenPVMatched_Pt2"]->Fill(recopvs[i].ptsq);
0824     } else {
0825       double pt2 = recopvs[i].ptsq;
0826       mes_[label]["RecoAssoc2GenPVNotMatched_Pt2"]->Fill(pt2);
0827       // Subtract hard-scatter track pt2 from the pileup pt2
0828       double pt2_pu = pt2 - vtx_sumpt2_sigmatched[i];
0829       mes_[label]["RecoAssoc2GenPVNotMatched_GenPVTracksRemoved_Pt2"]->Fill(pt2_pu);
0830     }
0831   }
0832   if (!signal_is_highest_pt && genpv_position_in_reco_collection >= 0)
0833     mes_[label]["RecoAssoc2GenPVMatchedNotHighest_Pt2"]->Fill(recopvs[genpv_position_in_reco_collection].ptsq);
0834 }
0835 
0836 /* Extract information form TrackingParticles/TrackingVertex and fill
0837  * the helper class simPrimaryVertex with proper generation-level
0838  * information */
0839 std::vector<PrimaryVertexAnalyzer4PUSlimmed::simPrimaryVertex> PrimaryVertexAnalyzer4PUSlimmed::getSimPVs(
0840     const edm::Handle<TrackingVertexCollection>& tVC) {
0841   std::vector<PrimaryVertexAnalyzer4PUSlimmed::simPrimaryVertex> simpv;
0842   int current_event = -1;
0843 
0844   if (verbose_) {
0845     std::cout << "getSimPVs TrackingVertexCollection " << std::endl;
0846   }
0847 
0848   for (TrackingVertexCollection::const_iterator v = tVC->begin(); v != tVC->end(); ++v) {
0849     if (verbose_) {
0850       std::cout << "BunchX.EventId: " << v->eventId().bunchCrossing() << "." << (v->eventId()).event()
0851                 << " Position: " << v->position() << " G4/HepMC Vertices: " << v->g4Vertices().size() << "/"
0852                 << v->genVertices().size() << "   t = " << v->position().t() * 1.e12
0853                 << "    == 0:" << (v->position().t() > 0) << std::endl;
0854       for (TrackingVertex::g4v_iterator gv = v->g4Vertices_begin(); gv != v->g4Vertices_end(); gv++) {
0855         std::cout << *gv << std::endl;
0856       }
0857       std::cout << "----------" << std::endl;
0858 
0859     }  // end of verbose_ session
0860 
0861     // I'd rather change this and select only vertices that come from
0862     // BX=0.  We should keep only the first vertex from all the events
0863     // at BX=0.
0864     if (v->eventId().bunchCrossing() != 0)
0865       continue;
0866     if (v->eventId().event() != current_event) {
0867       current_event = v->eventId().event();
0868     } else {
0869       continue;
0870     }
0871     // TODO(rovere) is this really necessary?
0872     if (fabs(v->position().z()) > 1000)
0873       continue;  // skip funny junk vertices
0874 
0875     // could be a new vertex, check  all primaries found so far to avoid
0876     // multiple entries
0877     simPrimaryVertex sv(v->position().x(), v->position().y(), v->position().z());
0878     sv.eventId = v->eventId();
0879     sv.sim_vertex = TrackingVertexRef(tVC, std::distance(tVC->begin(), v));
0880 
0881     for (TrackingParticleRefVector::iterator iTrack = v->daughterTracks_begin(); iTrack != v->daughterTracks_end();
0882          ++iTrack) {
0883       // TODO(rovere) isn't it always the case? Is it really worth
0884       // checking this out?
0885       // sv.eventId = (**iTrack).eventId();
0886       assert((**iTrack).eventId().bunchCrossing() == 0);
0887     }
0888     // TODO(rovere) maybe get rid of this old logic completely ... ?
0889     simPrimaryVertex* vp = nullptr;  // will become non-NULL if a vertex
0890                                      // is found and then point to it
0891     for (std::vector<simPrimaryVertex>::iterator v0 = simpv.begin(); v0 != simpv.end(); v0++) {
0892       if ((sv.eventId == v0->eventId) && (fabs(sv.x - v0->x) < 1e-5) && (fabs(sv.y - v0->y) < 1e-5) &&
0893           (fabs(sv.z - v0->z) < 1e-5)) {
0894         vp = &(*v0);
0895         break;
0896       }
0897     }
0898     if (!vp) {
0899       // this is a new vertex, add it to the list of sim-vertices
0900       simpv.push_back(sv);
0901       vp = &simpv.back();
0902       if (verbose_) {
0903         std::cout << "this is a new vertex " << sv.eventId.event() << "   " << sv.x << " " << sv.y << " " << sv.z
0904                   << std::endl;
0905       }
0906     } else {
0907       if (verbose_) {
0908         std::cout << "this is not a new vertex" << sv.x << " " << sv.y << " " << sv.z << std::endl;
0909       }
0910     }
0911 
0912     // Loop over daughter track(s) as Tracking Particles
0913     for (TrackingVertex::tp_iterator iTP = v->daughterTracks_begin(); iTP != v->daughterTracks_end(); ++iTP) {
0914       auto momentum = (*(*iTP)).momentum();
0915       const reco::Track* matched_best_reco_track = nullptr;
0916       double match_quality = -1;
0917       if (use_only_charged_tracks_ && (**iTP).charge() == 0)
0918         continue;
0919       if (s2r_->find(*iTP) != s2r_->end()) {
0920         matched_best_reco_track = (*s2r_)[*iTP][0].first.get();
0921         match_quality = (*s2r_)[*iTP][0].second;
0922       }
0923       if (verbose_) {
0924         std::cout << "  Daughter momentum:      " << momentum;
0925         std::cout << "  Daughter type     " << (*(*iTP)).pdgId();
0926         std::cout << "  matched: " << (matched_best_reco_track != nullptr);
0927         std::cout << "  match-quality: " << match_quality;
0928         std::cout << std::endl;
0929       }
0930       vp->ptot.setPx(vp->ptot.x() + momentum.x());
0931       vp->ptot.setPy(vp->ptot.y() + momentum.y());
0932       vp->ptot.setPz(vp->ptot.z() + momentum.z());
0933       vp->ptot.setE(vp->ptot.e() + (**iTP).energy());
0934       vp->ptsq += ((**iTP).pt() * (**iTP).pt());
0935       // TODO(rovere) only select charged sim-particles? If so, maybe
0936       // put it as a configuration parameter?
0937       if (matched_best_reco_track) {
0938         vp->num_matched_reco_tracks++;
0939         vp->average_match_quality += match_quality;
0940       }
0941       // TODO(rovere) get rid of cuts on sim-tracks
0942       // TODO(rovere) be consistent between simulated tracks and
0943       // reconstructed tracks selection
0944       // count relevant particles
0945       if (((**iTP).pt() > 0.2) && (fabs((**iTP).eta()) < 2.5) && (**iTP).charge() != 0) {
0946         vp->nGenTrk++;
0947       }
0948     }  // End of for loop on daughters sim-particles
0949     if (vp->num_matched_reco_tracks)
0950       vp->average_match_quality /= static_cast<float>(vp->num_matched_reco_tracks);
0951     if (verbose_) {
0952       std::cout << "average number of associated tracks: "
0953                 << vp->num_matched_reco_tracks / static_cast<float>(vp->nGenTrk)
0954                 << " with average quality: " << vp->average_match_quality << std::endl;
0955     }
0956   }  // End of for loop on tracking vertices
0957 
0958   if (verbose_) {
0959     std::cout << "------- PrimaryVertexAnalyzer4PUSlimmed simPVs from "
0960                  "TrackingVertices "
0961                  "-------"
0962               << std::endl;
0963     for (std::vector<simPrimaryVertex>::iterator v0 = simpv.begin(); v0 != simpv.end(); v0++) {
0964       std::cout << "z=" << v0->z << "  event=" << v0->eventId.event() << std::endl;
0965     }
0966     std::cout << "-----------------------------------------------" << std::endl;
0967   }  // End of for summary on discovered simulated primary vertices.
0968 
0969   // In case of no simulated vertices, break here
0970   if (simpv.empty())
0971     return simpv;
0972 
0973   // Now compute the closest distance in z between all simulated vertex
0974   // first initialize
0975   auto prev_z = simpv.back().z;
0976   for (simPrimaryVertex& vsim : simpv) {
0977     vsim.closest_vertex_distance_z = std::abs(vsim.z - prev_z);
0978     prev_z = vsim.z;
0979   }
0980   // then calculate
0981   for (std::vector<simPrimaryVertex>::iterator vsim = simpv.begin(); vsim != simpv.end(); vsim++) {
0982     std::vector<simPrimaryVertex>::iterator vsim2 = vsim;
0983     vsim2++;
0984     for (; vsim2 != simpv.end(); vsim2++) {
0985       double distance = std::abs(vsim->z - vsim2->z);
0986       // need both to be complete
0987       vsim->closest_vertex_distance_z = std::min(vsim->closest_vertex_distance_z, distance);
0988       vsim2->closest_vertex_distance_z = std::min(vsim2->closest_vertex_distance_z, distance);
0989     }
0990   }
0991   return simpv;
0992 }
0993 
0994 /* Extract information form recoVertex and fill the helper class
0995  * recoPrimaryVertex with proper reco-level information */
0996 std::vector<PrimaryVertexAnalyzer4PUSlimmed::recoPrimaryVertex> PrimaryVertexAnalyzer4PUSlimmed::getRecoPVs(
0997     const edm::Handle<edm::View<reco::Vertex>>& tVC) {
0998   std::vector<PrimaryVertexAnalyzer4PUSlimmed::recoPrimaryVertex> recopv;
0999 
1000   if (verbose_) {
1001     std::cout << "getRecoPVs TrackingVertexCollection " << std::endl;
1002   }
1003 
1004   for (auto v = tVC->begin(); v != tVC->end(); ++v) {
1005     if (verbose_) {
1006       std::cout << " Position: " << v->position() << std::endl;
1007     }
1008 
1009     // Skip junk vertices
1010     if (fabs(v->z()) > 1000)
1011       continue;
1012     if (v->isFake() || !v->isValid())
1013       continue;
1014 
1015     recoPrimaryVertex sv(v->position().x(), v->position().y(), v->position().z());
1016     sv.recVtx = &(*v);
1017     sv.recVtxRef = reco::VertexBaseRef(tVC, std::distance(tVC->begin(), v));
1018     // this is a new vertex, add it to the list of reco-vertices
1019     recopv.push_back(sv);
1020     PrimaryVertexAnalyzer4PUSlimmed::recoPrimaryVertex* vp = &recopv.back();
1021 
1022     // Loop over daughter track(s)
1023     for (auto iTrack = v->tracks_begin(); iTrack != v->tracks_end(); ++iTrack) {
1024       auto momentum = (*(*iTrack)).innerMomentum();
1025       // TODO(rovere) better handle the pixelVertices, whose tracks
1026       // do not have the innerMomentum defined. This is a temporary
1027       // hack to overcome this problem.
1028       if (momentum.mag2() == 0)
1029         momentum = (*(*iTrack)).momentum();
1030       if (verbose_) {
1031         std::cout << "  Daughter momentum:      " << momentum;
1032         std::cout << std::endl;
1033       }
1034       vp->pt += std::sqrt(momentum.perp2());
1035       vp->ptsq += (momentum.perp2());
1036       vp->nRecoTrk++;
1037 
1038       auto matched = r2s_->find(*iTrack);
1039       if (matched != r2s_->end()) {
1040         vp->num_matched_sim_tracks++;
1041       }
1042 
1043     }  // End of for loop on daughters reconstructed tracks
1044   }    // End of for loop on tracking vertices
1045 
1046   if (verbose_) {
1047     std::cout << "------- PrimaryVertexAnalyzer4PUSlimmed recoPVs from "
1048                  "VertexCollection "
1049                  "-------"
1050               << std::endl;
1051     for (std::vector<recoPrimaryVertex>::iterator v0 = recopv.begin(); v0 != recopv.end(); v0++) {
1052       std::cout << "z=" << v0->z << std::endl;
1053     }
1054     std::cout << "-----------------------------------------------" << std::endl;
1055   }  // End of for summary on reconstructed primary vertices.
1056 
1057   // In case of no reco vertices, break here
1058   if (recopv.empty())
1059     return recopv;
1060 
1061   // Now compute the closest distance in z between all reconstructed vertex
1062   // first initialize
1063   auto prev_z = recopv.back().z;
1064   for (recoPrimaryVertex& vreco : recopv) {
1065     vreco.closest_vertex_distance_z = std::abs(vreco.z - prev_z);
1066     prev_z = vreco.z;
1067   }
1068   for (std::vector<recoPrimaryVertex>::iterator vreco = recopv.begin(); vreco != recopv.end(); vreco++) {
1069     std::vector<recoPrimaryVertex>::iterator vreco2 = vreco;
1070     vreco2++;
1071     for (; vreco2 != recopv.end(); vreco2++) {
1072       double distance = std::abs(vreco->z - vreco2->z);
1073       // need both to be complete
1074       vreco->closest_vertex_distance_z = std::min(vreco->closest_vertex_distance_z, distance);
1075       vreco2->closest_vertex_distance_z = std::min(vreco2->closest_vertex_distance_z, distance);
1076     }
1077   }
1078   return recopv;
1079 }
1080 
1081 void PrimaryVertexAnalyzer4PUSlimmed::resetSimPVAssociation(std::vector<simPrimaryVertex>& simpv) {
1082   for (auto& v : simpv) {
1083     v.rec_vertices.clear();
1084   }
1085 }
1086 
1087 // ------------ method called to produce the data  ------------
1088 void PrimaryVertexAnalyzer4PUSlimmed::matchSim2RecoVertices(std::vector<simPrimaryVertex>& simpv,
1089                                                             const reco::VertexSimToRecoCollection& vertex_s2r) {
1090   if (verbose_) {
1091     std::cout << "PrimaryVertexAnalyzer4PUSlimmed::matchSim2RecoVertices " << std::endl;
1092   }
1093   for (std::vector<simPrimaryVertex>::iterator vsim = simpv.begin(); vsim != simpv.end(); vsim++) {
1094     auto matched = vertex_s2r.find(vsim->sim_vertex);
1095     if (matched != vertex_s2r.end()) {
1096       for (const auto& vertexRefQuality : matched->val) {
1097         vsim->rec_vertices.push_back(&(*(vertexRefQuality.first)));
1098       }
1099     }
1100 
1101     if (verbose_) {
1102       if (!vsim->rec_vertices.empty()) {
1103         for (auto const& v : vsim->rec_vertices) {
1104           std::cout << "Found a matching vertex for genVtx " << vsim->z << " at " << v->z()
1105                     << " with sign: " << fabs(v->z() - vsim->z) / v->zError() << std::endl;
1106         }
1107       } else {
1108         std::cout << "No matching vertex for " << vsim->z << std::endl;
1109       }
1110     }
1111   }  // end for loop on simulated vertices
1112   if (verbose_) {
1113     std::cout << "Done with matching sim vertices" << std::endl;
1114   }
1115 }
1116 
1117 void PrimaryVertexAnalyzer4PUSlimmed::matchReco2SimVertices(std::vector<recoPrimaryVertex>& recopv,
1118                                                             const reco::VertexRecoToSimCollection& vertex_r2s,
1119                                                             const std::vector<simPrimaryVertex>& simpv) {
1120   for (std::vector<recoPrimaryVertex>::iterator vrec = recopv.begin(); vrec != recopv.end(); vrec++) {
1121     auto matched = vertex_r2s.find(vrec->recVtxRef);
1122     if (matched != vertex_r2s.end()) {
1123       for (const auto& vertexRefQuality : matched->val) {
1124         const auto tvPtr = &(*(vertexRefQuality.first));
1125         vrec->sim_vertices.push_back(tvPtr);
1126       }
1127 
1128       for (const TrackingVertex* tv : vrec->sim_vertices) {
1129         // Set pointers to internal simVertex objects
1130         for (const auto& vv : simpv) {
1131           if (&(*(vv.sim_vertex)) == tv) {
1132             vrec->sim_vertices_internal.push_back(&vv);
1133             continue;
1134           }
1135         }
1136 
1137         // Calculate number of shared tracks
1138         vrec->sim_vertices_num_shared_tracks.push_back(calculateVertexSharedTracks(*(vrec->recVtx), *tv, *r2s_));
1139       }
1140     }
1141 
1142     if (verbose_) {
1143       for (auto v : vrec->sim_vertices) {
1144         std::cout << "Found a matching vertex for reco: " << vrec->z << " at gen:" << v->position().z()
1145                   << " with sign: " << fabs(vrec->z - v->position().z()) / vrec->recVtx->zError() << std::endl;
1146       }
1147     }
1148   }  // end for loop on reconstructed vertices
1149 }
1150 
1151 void PrimaryVertexAnalyzer4PUSlimmed::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
1152   using edm::Handle;
1153   using edm::View;
1154   using std::cout;
1155   using std::endl;
1156   using std::vector;
1157   using namespace reco;
1158 
1159   std::vector<float> pileUpInfo_z;
1160 
1161   // get the pileup information
1162   edm::Handle<std::vector<PileupSummaryInfo>> puinfoH;
1163   if (iEvent.getByToken(vecPileupSummaryInfoToken_, puinfoH)) {
1164     for (auto const& pu_info : *puinfoH.product()) {
1165       if (do_generic_sim_plots_) {
1166         mes_["root_folder"]["GenVtx_vs_BX"]->Fill(pu_info.getBunchCrossing(), pu_info.getPU_NumInteractions());
1167       }
1168       if (pu_info.getBunchCrossing() == 0) {
1169         pileUpInfo_z = pu_info.getPU_zpositions();
1170         if (verbose_) {
1171           for (auto const& p : pileUpInfo_z) {
1172             std::cout << "PileUpInfo on Z vertex: " << p << std::endl;
1173           }
1174         }
1175         break;
1176       }
1177     }
1178   }
1179 
1180   edm::Handle<TrackingParticleCollection> TPCollectionH;
1181   iEvent.getByToken(trackingParticleCollectionToken_, TPCollectionH);
1182   if (!TPCollectionH.isValid())
1183     edm::LogWarning("PrimaryVertexAnalyzer4PUSlimmed") << "TPCollectionH is not valid";
1184 
1185   edm::Handle<TrackingVertexCollection> TVCollectionH;
1186   iEvent.getByToken(trackingVertexCollectionToken_, TVCollectionH);
1187   if (!TVCollectionH.isValid())
1188     edm::LogWarning("PrimaryVertexAnalyzer4PUSlimmed") << "TVCollectionH is not valid";
1189 
1190   // TODO(rovere) the idea is to put in case a track-selector in front
1191   // of this module and then use its label to get the selected tracks
1192   // out of the event instead of making an hard-coded selection in the
1193   // code.
1194 
1195   edm::Handle<reco::SimToRecoCollection> simToRecoH;
1196   iEvent.getByToken(simToRecoAssociationToken_, simToRecoH);
1197   if (simToRecoH.isValid())
1198     s2r_ = simToRecoH.product();
1199   else
1200     edm::LogWarning("PrimaryVertexAnalyzer4PUSlimmed") << "simToRecoH is not valid";
1201 
1202   edm::Handle<reco::RecoToSimCollection> recoToSimH;
1203   iEvent.getByToken(recoToSimAssociationToken_, recoToSimH);
1204   if (recoToSimH.isValid())
1205     r2s_ = recoToSimH.product();
1206   else
1207     edm::LogWarning("PrimaryVertexAnalyzer4PUSlimmed") << "recoToSimH is not valid";
1208 
1209   // Vertex associator
1210   edm::Handle<reco::VertexToTrackingVertexAssociator> vertexAssociatorH;
1211   iEvent.getByToken(vertexAssociatorToken_, vertexAssociatorH);
1212   if (!vertexAssociatorH.isValid()) {
1213     edm::LogWarning("PrimaryVertexAnalyzer4PUSlimmed") << "vertexAssociatorH is not valid";
1214     return;
1215   }
1216   const reco::VertexToTrackingVertexAssociator& vertexAssociator = *(vertexAssociatorH.product());
1217 
1218   std::vector<simPrimaryVertex> simpv;  // a list of simulated primary
1219                                         // MC vertices
1220   // TODO(rovere) use move semantic?
1221   simpv = getSimPVs(TVCollectionH);
1222   // TODO(rovere) 1 vertex is not, by definition, pileup, and should
1223   // probably be subtracted?
1224   int kind_of_signal_vertex = 0;
1225   int num_pileup_vertices = simpv.size();
1226   if (do_generic_sim_plots_)
1227     mes_["root_folder"]["GenAllV_NumVertices"]->Fill(simpv.size());
1228   bool signal_is_highest_pt =
1229       std::max_element(simpv.begin(), simpv.end(), [](const simPrimaryVertex& lhs, const simPrimaryVertex& rhs) {
1230         return lhs.ptsq < rhs.ptsq;
1231       }) == simpv.begin();
1232   kind_of_signal_vertex |= (signal_is_highest_pt << HIGHEST_PT);
1233   if (do_generic_sim_plots_) {
1234     mes_["root_folder"]["SignalIsHighestPt2"]->Fill(signal_is_highest_pt ? 1. : 0.);
1235     computePairDistance(simpv, mes_["root_folder"]["GenAllV_PairDistanceZ"]);
1236   }
1237 
1238   int label_index = -1;
1239   for (size_t iToken = 0, endToken = reco_vertex_collection_tokens_.size(); iToken < endToken; ++iToken) {
1240     auto const& vertex_token = reco_vertex_collection_tokens_[iToken];
1241     std::vector<recoPrimaryVertex> recopv;  // a list of reconstructed
1242                                             // primary MC vertices
1243     std::string label = reco_vertex_collections_[++label_index].label();
1244     edm::Handle<edm::View<reco::Vertex>> recVtxs;
1245     if (!iEvent.getByToken(vertex_token, recVtxs)) {
1246       if (!errorPrintedForColl_[iToken]) {
1247         edm::LogWarning("PrimaryVertexAnalyzer4PUSlimmed")
1248             << "Skipping vertex collection: " << label << " since it is missing.";
1249         errorPrintedForColl_[iToken] = true;
1250       }
1251       continue;
1252     }
1253 
1254     {
1255       // check upfront that refs to track are (likely) to be valid
1256       bool ok = true;
1257       for (const auto& v : *recVtxs) {
1258         if (v.tracksSize() > 0) {
1259           const auto& ref = v.trackRefAt(0);
1260           if (ref.isNull() || !ref.isAvailable()) {
1261             if (!errorPrintedForColl_[iToken]) {
1262               edm::LogWarning("PrimaryVertexAnalyzer4PUSlimmed")
1263                   << "Skipping vertex collection: " << label
1264                   << " since likely the track collection the vertex has refs pointing to is missing (at least the "
1265                      "first TrackBaseRef is null or not available)";
1266               errorPrintedForColl_[iToken] = true;
1267             }
1268             ok = false;
1269           }
1270         }
1271       }
1272       if (!ok)
1273         continue;
1274     }
1275 
1276     reco::VertexRecoToSimCollection vertex_r2s = vertexAssociator.associateRecoToSim(recVtxs, TVCollectionH);
1277     reco::VertexSimToRecoCollection vertex_s2r = vertexAssociator.associateSimToReco(recVtxs, TVCollectionH);
1278 
1279     resetSimPVAssociation(simpv);
1280     matchSim2RecoVertices(simpv, vertex_s2r);
1281     recopv = getRecoPVs(recVtxs);
1282     computePairDistance(recopv, mes_[label]["RecoAllAssoc2Gen_PairDistanceZ"]);
1283     matchReco2SimVertices(recopv, vertex_r2s, simpv);
1284 
1285     int num_total_gen_vertices_assoc2reco = 0;
1286     int num_total_reco_vertices_assoc2gen = 0;
1287     int num_total_gen_vertices_multiassoc2reco = 0;
1288     int num_total_reco_vertices_multiassoc2gen = 0;
1289     int num_total_reco_vertices_duplicate = 0;
1290     int genpv_position_in_reco_collection = -1;
1291     for (auto const& v : simpv) {
1292       float mistag = 1.;
1293       // TODO(rovere) put selectors here in front of fill* methods.
1294       if (v.eventId.event() == 0) {
1295         if (!recVtxs->empty() && std::find(v.rec_vertices.begin(), v.rec_vertices.end(), &((*recVtxs.product())[0])) !=
1296                                      v.rec_vertices.end()) {
1297           mistag = 0.;
1298           kind_of_signal_vertex |= (1 << IS_ASSOC2FIRST_RECO);
1299         } else {
1300           if (!v.rec_vertices.empty()) {
1301             kind_of_signal_vertex |= (1 << IS_ASSOC2ANY_RECO);
1302           }
1303         }
1304         mes_[label]["KindOfSignalPV"]->Fill(kind_of_signal_vertex);
1305         mes_[label]["MisTagRate"]->Fill(mistag);
1306         mes_[label]["MisTagRate_vs_PU"]->Fill(simpv.size(), mistag);
1307         mes_[label]["MisTagRate_vs_sum-pt2"]->Fill(v.ptsq, mistag);
1308         mes_[label]["MisTagRate_vs_Z"]->Fill(v.z, mistag);
1309         mes_[label]["MisTagRate_vs_R"]->Fill(v.r, mistag);
1310         mes_[label]["MisTagRate_vs_NumTracks"]->Fill(v.nGenTrk, mistag);
1311         if (signal_is_highest_pt) {
1312           mes_[label]["MisTagRateSignalIsHighest"]->Fill(mistag);
1313           mes_[label]["MisTagRateSignalIsHighest_vs_PU"]->Fill(simpv.size(), mistag);
1314           mes_[label]["MisTagRateSignalIsHighest_vs_sum-pt2"]->Fill(v.ptsq, mistag);
1315           mes_[label]["MisTagRateSignalIsHighest_vs_Z"]->Fill(v.z, mistag);
1316           mes_[label]["MisTagRateSignalIsHighest_vs_R"]->Fill(v.r, mistag);
1317           mes_[label]["MisTagRateSignalIsHighest_vs_NumTracks"]->Fill(v.nGenTrk, mistag);
1318         } else {
1319           mes_[label]["MisTagRateSignalIsNotHighest"]->Fill(mistag);
1320           mes_[label]["MisTagRateSignalIsNotHighest_vs_PU"]->Fill(simpv.size(), mistag);
1321           mes_[label]["MisTagRateSignalIsNotHighest_vs_sum-pt2"]->Fill(v.ptsq, mistag);
1322           mes_[label]["MisTagRateSignalIsNotHighest_vs_Z"]->Fill(v.z, mistag);
1323           mes_[label]["MisTagRateSignalIsNotHighest_vs_R"]->Fill(v.r, mistag);
1324           mes_[label]["MisTagRateSignalIsNotHighest_vs_NumTracks"]->Fill(v.nGenTrk, mistag);
1325         }
1326         // Now check at which location the Simulated PV has been
1327         // reconstructed in the primary vertex collection
1328         // at-hand. Mark it with fake index -1 if it was not
1329         // reconstructed at all.
1330 
1331         auto iv = (*recVtxs.product()).begin();
1332         for (int pv_position_in_reco_collection = 0; iv != (*recVtxs.product()).end();
1333              ++pv_position_in_reco_collection, ++iv) {
1334           if (std::find(v.rec_vertices.begin(), v.rec_vertices.end(), &(*iv)) != v.rec_vertices.end()) {
1335             mes_[label]["TruePVLocationIndex"]->Fill(pv_position_in_reco_collection);
1336             const bool genPVMatchedToRecoPV = (pv_position_in_reco_collection == 0);
1337             mes_[label]["TruePVLocationIndexCumulative"]->Fill(genPVMatchedToRecoPV ? 0 : 1);
1338 
1339             if (signal_is_highest_pt) {
1340               mes_[label]["TruePVLocationIndexSignalIsHighest"]->Fill(pv_position_in_reco_collection);
1341             } else {
1342               mes_[label]["TruePVLocationIndexSignalIsNotHighest"]->Fill(pv_position_in_reco_collection);
1343             }
1344 
1345             fillRecoAssociatedGenPVHistograms(label, v, genPVMatchedToRecoPV);
1346             if (genPVMatchedToRecoPV) {
1347               auto pv = recopv[0];
1348               assert(pv.recVtx == &(*iv));
1349               fillResolutionAndPullHistograms(label, num_pileup_vertices, pv, true);
1350             }
1351             genpv_position_in_reco_collection = pv_position_in_reco_collection;
1352             break;
1353           }
1354         }
1355 
1356         // If we reached the end, it means that the Simulated PV has not
1357         // been associated to any reconstructed vertex: mark it as
1358         // missing in the reconstructed vertex collection using the fake
1359         // index -1.
1360         if (iv == (*recVtxs.product()).end()) {
1361           mes_[label]["TruePVLocationIndex"]->Fill(-1.);
1362           mes_[label]["TruePVLocationIndexCumulative"]->Fill(-1.);
1363           if (signal_is_highest_pt)
1364             mes_[label]["TruePVLocationIndexSignalIsHighest"]->Fill(-1.);
1365           else
1366             mes_[label]["TruePVLocationIndexSignalIsNotHighest"]->Fill(-1.);
1367         }
1368       }
1369 
1370       if (!v.rec_vertices.empty())
1371         num_total_gen_vertices_assoc2reco++;
1372       if (v.rec_vertices.size() > 1)
1373         num_total_gen_vertices_multiassoc2reco++;
1374       // No need to N-tplicate the Gen-related cumulative histograms:
1375       // fill them only at the first iteration
1376       if (do_generic_sim_plots_ && label_index == 0)
1377         fillGenericGenVertexHistograms(v);
1378       fillRecoAssociatedGenVertexHistograms(label, v);
1379     }
1380     calculatePurityAndFillHistograms(label, recopv, genpv_position_in_reco_collection, signal_is_highest_pt);
1381 
1382     mes_[label]["GenAllAssoc2Reco_NumVertices"]->Fill(simpv.size(), simpv.size());
1383     mes_[label]["GenAllAssoc2RecoMatched_NumVertices"]->Fill(simpv.size(), num_total_gen_vertices_assoc2reco);
1384     mes_[label]["GenAllAssoc2RecoMultiMatched_NumVertices"]->Fill(simpv.size(), num_total_gen_vertices_multiassoc2reco);
1385     for (auto& v : recopv) {
1386       fillGenAssociatedRecoVertexHistograms(label, num_pileup_vertices, v);
1387       if (!v.sim_vertices.empty()) {
1388         num_total_reco_vertices_assoc2gen++;
1389         if (v.sim_vertices_internal[0]->rec_vertices.size() > 1) {
1390           num_total_reco_vertices_duplicate++;
1391         }
1392       }
1393       if (v.sim_vertices.size() > 1)
1394         num_total_reco_vertices_multiassoc2gen++;
1395     }
1396     mes_[label]["RecoAllAssoc2Gen_NumVertices"]->Fill(recopv.size(), recopv.size());
1397     mes_[label]["RecoAllAssoc2GenMatched_NumVertices"]->Fill(recopv.size(), num_total_reco_vertices_assoc2gen);
1398     mes_[label]["RecoAllAssoc2GenMultiMatched_NumVertices"]->Fill(recopv.size(),
1399                                                                   num_total_reco_vertices_multiassoc2gen);
1400     mes_[label]["RecoAllAssoc2MultiMatchedGen_NumVertices"]->Fill(recopv.size(), num_total_reco_vertices_duplicate);
1401     mes_[label]["RecoVtx_vs_GenVtx"]->Fill(simpv.size(), recopv.size());
1402     mes_[label]["MatchedRecoVtx_vs_GenVtx"]->Fill(simpv.size(), num_total_reco_vertices_assoc2gen);
1403   }
1404 }  // end of analyze
1405 
1406 template <class T>
1407 void PrimaryVertexAnalyzer4PUSlimmed::computePairDistance(const T& collection, MonitorElement* me) {
1408   for (unsigned int i = 0; i < collection.size(); ++i) {
1409     for (unsigned int j = i + 1; j < collection.size(); ++j) {
1410       me->Fill(std::abs(collection[i].z - collection[j].z));
1411     }
1412   }
1413 }
[ 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!