File indexing completed on 2025-05-20 22:24:40
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017 #include "TMath.h"
0018 #include "TFile.h"
0019 #include "TH1I.h"
0020 #include "TH1D.h"
0021 #include "TH2D.h"
0022 #include "TProfile.h"
0023 #include "TProfile2D.h"
0024
0025
0026 #include <memory>
0027 #include <fmt/format.h>
0028 #include <fmt/printf.h>
0029
0030
0031 #include "CommonTools/Statistics/interface/ChiSquaredProbability.h"
0032 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0033 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0034 #include "DataFormats/Common/interface/Association.h"
0035 #include "DataFormats/Common/interface/RefToBase.h"
0036 #include "DataFormats/Common/interface/ValueMap.h"
0037 #include "DataFormats/TrackReco/interface/Track.h"
0038 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0039 #include "DataFormats/VertexReco/interface/Vertex.h"
0040 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0041 #include "FWCore/Framework/interface/Event.h"
0042 #include "FWCore/Framework/interface/Frameworkfwd.h"
0043 #include "FWCore/Framework/interface/MakerMacros.h"
0044 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0045 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0046 #include "FWCore/ServiceRegistry/interface/Service.h"
0047 #include "FWCore/Utilities/interface/InputTag.h"
0048 #include "FWCore/Utilities/interface/isFinite.h"
0049
0050
0051
0052
0053
0054 using reco::TrackCollection;
0055
0056 namespace gpVertexAnalyzer {
0057 void setBinLog(TAxis *axis) {
0058 int bins = axis->GetNbins();
0059 float from = axis->GetXmin();
0060 float to = axis->GetXmax();
0061 float width = (to - from) / bins;
0062 std::vector<float> new_bins(bins + 1, 0);
0063 for (int i = 0; i <= bins; i++) {
0064 new_bins[i] = TMath::Power(10, from + i * width);
0065 }
0066 axis->Set(bins, new_bins.data());
0067 }
0068
0069 void setBinLogX(TH1 *h) {
0070 TAxis *axis = h->GetXaxis();
0071 setBinLog(axis);
0072 }
0073
0074 void setBinLogY(TH1 *h) {
0075 TAxis *axis = h->GetYaxis();
0076 setBinLog(axis);
0077 }
0078
0079 template <typename... Args>
0080 TProfile *makeProfileIfLog(const edm::Service<TFileService> &fs, bool logx, bool logy, Args &&...args) {
0081 auto prof = fs->make<TProfile>(std::forward<Args>(args)...);
0082 if (logx)
0083 setBinLogX(prof);
0084 if (logy)
0085 setBinLogY(prof);
0086 return prof;
0087 }
0088
0089 template <typename... Args>
0090 TH1D *makeTH1IfLog(const edm::Service<TFileService> &fs, bool logx, bool logy, Args &&...args) {
0091 auto h1 = fs->make<TH1D>(std::forward<Args>(args)...);
0092 if (logx)
0093 setBinLogX(h1);
0094 if (logy)
0095 setBinLogY(h1);
0096 return h1;
0097 }
0098 }
0099
0100 class GeneralPurposeVertexAnalyzer : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0101 public:
0102 explicit GeneralPurposeVertexAnalyzer(const edm::ParameterSet &);
0103 ~GeneralPurposeVertexAnalyzer() override = default;
0104
0105 static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0106
0107 private:
0108 void pvTracksPlots(const reco::Vertex &v);
0109 void vertexPlots(const reco::Vertex &v, const reco::BeamSpot &beamSpot, int i);
0110 template <typename T, typename... Args>
0111 T *book(const Args &...args) const;
0112 void beginJob() override;
0113 void analyze(const edm::Event &, const edm::EventSetup &) override;
0114
0115
0116 edm::Service<TFileService> fs_;
0117
0118 struct IPMonitoring {
0119 std::string varname_;
0120 float pTcut_;
0121 TH1D *IP_, *IPErr_, *IPPull_;
0122 TProfile *IPVsPhi_, *IPVsEta_, *IPVsPt_;
0123 TProfile *IPErrVsPhi_, *IPErrVsEta_, *IPErrVsPt_;
0124 TProfile2D *IPVsEtaVsPhi_, *IPErrVsEtaVsPhi_;
0125
0126 void bookIPMonitor(const edm::ParameterSet &, const edm::Service<TFileService> fs);
0127 void fillIPMonitor(const edm::RefToBase<reco::Track> &tk, const math::XYZPoint &pv);
0128 };
0129
0130 const edm::ParameterSet conf_;
0131 const int ndof_;
0132 bool errorPrinted_;
0133
0134 const edm::InputTag vertexInputTag_, beamSpotInputTag_;
0135 const edm::EDGetTokenT<reco::VertexCollection> vertexToken_;
0136 using VertexScore = edm::ValueMap<float>;
0137 const edm::EDGetTokenT<VertexScore> scoreToken_;
0138 const edm::EDGetTokenT<reco::BeamSpot> beamspotToken_;
0139
0140 static constexpr int cmToUm = 10000;
0141
0142 const double vposx_;
0143 const double vposy_;
0144 const int tkNoBin_;
0145 const double tkNoMin_;
0146 const double tkNoMax_;
0147
0148 const int dxyBin_;
0149 const double dxyMin_;
0150 const double dxyMax_;
0151
0152 const int dzBin_;
0153 const double dzMin_;
0154 const double dzMax_;
0155
0156 const int phiBin_;
0157 const int phiBin2D_;
0158 const double phiMin_;
0159 const double phiMax_;
0160
0161 const int etaBin_;
0162 const int etaBin2D_;
0163 const double etaMin_;
0164 const double etaMax_;
0165
0166
0167 TH1I *nbvtx, *nbgvtx;
0168 TH1D *nbtksinvtx[2], *trksWeight[2], *score[2];
0169 TH1D *tt[2];
0170 TH1D *xrec[2], *yrec[2], *zrec[2], *xDiff[2], *yDiff[2], *xerr[2], *yerr[2], *zerr[2];
0171 TH2D *xerrVsTrks[2], *yerrVsTrks[2], *zerrVsTrks[2];
0172 TH1D *ntracksVsZ[2];
0173 TH1D *vtxchi2[2], *vtxndf[2], *vtxprob[2], *nans[2];
0174 TH1D *type[2];
0175 TH1D *bsX, *bsY, *bsZ, *bsSigmaZ, *bsDxdz, *bsDydz, *bsBeamWidthX, *bsBeamWidthY, *bsType;
0176
0177 TH1D *trackpt, *sumpt, *ntracks, *weight, *chi2ndf, *chi2prob;
0178 TH1D *dxy2;
0179 TH1D *phi_pt1, *eta_pt1;
0180 TH1D *phi_pt10, *eta_pt10;
0181
0182
0183 IPMonitoring dxy_pt1;
0184 IPMonitoring dxy_pt10;
0185
0186 IPMonitoring dz_pt1;
0187 IPMonitoring dz_pt10;
0188 };
0189
0190 void GeneralPurposeVertexAnalyzer::IPMonitoring::bookIPMonitor(const edm::ParameterSet &config,
0191 const edm::Service<TFileService> fs) {
0192
0193 auto checkRange = [](const std::string &name, double min, double max) {
0194 if (min >= max) {
0195 throw cms::Exception("IPMonitoring")
0196 << "Invalid range for " << name << ": min (" << min << ") >= max (" << max << ")";
0197 }
0198 };
0199
0200 int VarBin = config.getParameter<int>(fmt::format("D{}Bin", varname_));
0201 double VarMin = config.getParameter<double>(fmt::format("D{}Min", varname_));
0202 double VarMax = config.getParameter<double>(fmt::format("D{}Max", varname_));
0203 checkRange(fmt::format("D{}", varname_), VarMin, VarMax);
0204
0205 int PhiBin = config.getParameter<int>("PhiBin");
0206 int PhiBin2D = config.getParameter<int>("PhiBin2D");
0207 double PhiMin = config.getParameter<double>("PhiMin");
0208 double PhiMax = config.getParameter<double>("PhiMax");
0209 checkRange("Phi", PhiMin, PhiMax);
0210
0211 int EtaBin = config.getParameter<int>("EtaBin");
0212 int EtaBin2D = config.getParameter<int>("EtaBin2D");
0213 double EtaMin = config.getParameter<double>("EtaMin");
0214 double EtaMax = config.getParameter<double>("EtaMax");
0215 checkRange("Eta", EtaMin, EtaMax);
0216
0217 int PtBin = config.getParameter<int>("PtBin");
0218 double PtMin = config.getParameter<double>("PtMin") * pTcut_;
0219 double PtMax = config.getParameter<double>("PtMax") * pTcut_;
0220 checkRange("Pt", PtMin, PtMax);
0221
0222 IP_ = fs->make<TH1D>(fmt::format("d{}_pt{}", varname_, pTcut_).c_str(),
0223 fmt::format("PV tracks (p_{{T}} > {} GeV) d_{{{}}} (#mum)", pTcut_, varname_).c_str(),
0224 VarBin,
0225 VarMin,
0226 VarMax);
0227
0228 IPErr_ = fs->make<TH1D>(fmt::format("d{}Err_pt{}", varname_, pTcut_).c_str(),
0229 fmt::format("PV tracks (p_{{T}} > {} GeV) d_{{{}}} error (#mum)", pTcut_, varname_).c_str(),
0230 100,
0231 0.,
0232 (varname_.find("xy") != std::string::npos) ? 2000. : 10000.);
0233
0234 IPPull_ = fs->make<TH1D>(
0235 fmt::format("d{}Pull_pt{}", varname_, pTcut_).c_str(),
0236 fmt::format("PV tracks (p_{{T}} > {} GeV) d_{{{}}}/#sigma_{{d_{{{}}}}}", pTcut_, varname_, varname_).c_str(),
0237 100,
0238 -5.,
0239 5.);
0240
0241 IPVsPhi_ =
0242 fs->make<TProfile>(fmt::format("d{}VsPhi_pt{}", varname_, pTcut_).c_str(),
0243 fmt::format("PV tracks (p_{{T}} > {}) d_{{{}}} VS track #phi", pTcut_, varname_).c_str(),
0244 PhiBin,
0245 PhiMin,
0246 PhiMax,
0247
0248 VarMin,
0249 VarMax);
0250 IPVsPhi_->SetXTitle("PV track (p_{T} > 1 GeV) #phi");
0251 IPVsPhi_->SetYTitle(fmt::format("PV tracks (p_{{T}} > {} GeV) d_{{{}}} (#mum)", pTcut_, varname_).c_str());
0252
0253 IPVsEta_ =
0254 fs->make<TProfile>(fmt::format("d{}VsEta_pt{}", varname_, pTcut_).c_str(),
0255 fmt::format("PV tracks (p_{{T}} > {}) d_{{{}}} VS track #eta", pTcut_, varname_).c_str(),
0256 EtaBin,
0257 EtaMin,
0258 EtaMax,
0259
0260 VarMin,
0261 VarMax);
0262 IPVsEta_->SetXTitle("PV track (p_{T} > 1 GeV) #eta");
0263 IPVsEta_->SetYTitle(fmt::format("PV tracks (p_{{T}} > {} GeV) d_{{{}}} (#mum)", pTcut_, varname_).c_str());
0264
0265 IPVsPt_ = gpVertexAnalyzer::makeProfileIfLog(
0266 fs,
0267 true,
0268 false,
0269 fmt::format("d{}VsPt_pt{}", varname_, pTcut_).c_str(),
0270 fmt::format("PV tracks (p_{{T}} > {}) d_{{{}}} VS track p_{{T}}", pTcut_, varname_).c_str(),
0271 PtBin,
0272 log10(PtMin),
0273 log10(PtMax),
0274 VarMin,
0275 VarMax,
0276 "");
0277 IPVsPt_->SetXTitle("PV track (p_{T} > 1 GeV) p_{T} [GeV]");
0278 IPVsPt_->SetYTitle(fmt::format("PV tracks (p_{{T}} > {} GeV) d_{{{}}} (#mum)", pTcut_, varname_).c_str());
0279
0280 IPErrVsPhi_ =
0281 fs->make<TProfile>(fmt::format("d{}ErrVsPhi_pt{}", varname_, pTcut_).c_str(),
0282 fmt::format("PV tracks (p_{{T}} > {}) d_{{{}}} error VS track #phi", pTcut_, varname_).c_str(),
0283 PhiBin,
0284 PhiMin,
0285 PhiMax,
0286
0287 0.,
0288 (varname_.find("xy") != std::string::npos) ? 100. : 200.);
0289 IPErrVsPhi_->SetXTitle("PV track (p_{T} > 1 GeV) #phi");
0290 IPErrVsPhi_->SetYTitle(fmt::format("PV tracks (p_{{T}} > {} GeV) d_{{{}}} error (#mum)", pTcut_, varname_).c_str());
0291
0292 IPErrVsEta_ =
0293 fs->make<TProfile>(fmt::format("d{}ErrVsEta_pt{}", varname_, pTcut_).c_str(),
0294 fmt::format("PV tracks (p_{{T}} > {}) d_{{{}}} error VS track #eta", pTcut_, varname_).c_str(),
0295 EtaBin,
0296 EtaMin,
0297 EtaMax,
0298
0299 0.,
0300 (varname_.find("xy") != std::string::npos) ? 100. : 200.);
0301 IPErrVsEta_->SetXTitle("PV track (p_{T} > 1 GeV) #eta");
0302 IPErrVsEta_->SetYTitle(fmt::format("PV tracks (p_{{T}} > {} GeV) d_{{{}}} error (#mum)", pTcut_, varname_).c_str());
0303
0304 IPErrVsPt_ = gpVertexAnalyzer::makeProfileIfLog(
0305 fs,
0306 true,
0307 false,
0308 fmt::format("d{}ErrVsPt_pt{}", varname_, pTcut_).c_str(),
0309 fmt::format("PV tracks (p_{{T}} > {}) d_{{{}}} error VS track p_{{T}}", pTcut_, varname_).c_str(),
0310 PtBin,
0311 log10(PtMin),
0312 log10(PtMax),
0313 VarMin,
0314 VarMax,
0315 "");
0316 IPErrVsPt_->SetXTitle("PV track (p_{T} > 1 GeV) p_{T} [GeV]");
0317 IPErrVsPt_->SetYTitle(fmt::format("PV tracks (p_{{T}} > {} GeV) d_{{{}}} error (#mum)", pTcut_, varname_).c_str());
0318
0319 IPVsEtaVsPhi_ = fs->make<TProfile2D>(
0320 fmt::format("d{}VsEtaVsPhi_pt{}", varname_, pTcut_).c_str(),
0321 fmt::format("PV tracks (p_{{T}} > {}) d_{{{}}} VS track #eta VS track #phi", pTcut_, varname_).c_str(),
0322 EtaBin2D,
0323 EtaMin,
0324 EtaMax,
0325 PhiBin2D,
0326 PhiMin,
0327 PhiMax,
0328
0329 VarMin,
0330 VarMax);
0331 IPVsEtaVsPhi_->SetXTitle("PV track (p_{T} > 1 GeV) #eta");
0332 IPVsEtaVsPhi_->SetYTitle("PV track (p_{T} > 1 GeV) #phi");
0333 IPVsEtaVsPhi_->SetZTitle(fmt::format("PV tracks (p_{{T}} > {} GeV) d_{{{}}} (#mum)", pTcut_, varname_).c_str());
0334
0335 IPErrVsEtaVsPhi_ = fs->make<TProfile2D>(
0336 fmt::format("d{}ErrVsEtaVsPhi_pt{}", varname_, pTcut_).c_str(),
0337 fmt::format("PV tracks (p_{{T}} > {}) d_{{{}}} error VS track #eta VS track #phi", pTcut_, varname_).c_str(),
0338 EtaBin2D,
0339 EtaMin,
0340 EtaMax,
0341 PhiBin2D,
0342 PhiMin,
0343 PhiMax,
0344
0345 0.,
0346 (varname_.find("xy") != std::string::npos) ? 100. : 200.);
0347 IPErrVsEtaVsPhi_->SetXTitle("PV track (p_{T} > 1 GeV) #eta");
0348 IPErrVsEtaVsPhi_->SetYTitle("PV track (p_{T} > 1 GeV) #phi");
0349 IPErrVsEtaVsPhi_->SetZTitle(
0350 fmt::format("PV tracks (p_{{T}} > {} GeV) d_{{{}}} error (#mum)", pTcut_, varname_).c_str());
0351 }
0352
0353 void GeneralPurposeVertexAnalyzer::IPMonitoring::fillIPMonitor(const edm::RefToBase<reco::Track> &tk,
0354 const math::XYZPoint &pv) {
0355 float value = 999.f, error = 9999.f;
0356
0357 if (varname_ == "xy") {
0358 value = tk->dxy(pv) * cmToUm;
0359 error = tk->dxyError() * cmToUm;
0360 } else if (varname_ == "z") {
0361 value = tk->dz(pv) * cmToUm;
0362 error = tk->dzError() * cmToUm;
0363 } else {
0364 throw cms::Exception("IPMonitoring") << "Unknown varname: " << varname_ << ". Expected 'xy' or 'z\'.";
0365 }
0366
0367 const float eta = tk->eta();
0368 const float phi = tk->phi();
0369 const float pt = tk->pt();
0370
0371 IP_->Fill(value);
0372 IPVsPhi_->Fill(phi, value);
0373 IPVsEta_->Fill(eta, value);
0374 IPVsPt_->Fill(pt, value);
0375 IPVsEtaVsPhi_->Fill(eta, phi, value);
0376
0377 IPErr_->Fill(error);
0378 IPPull_->Fill(value / error);
0379 IPErrVsPhi_->Fill(phi, error);
0380 IPErrVsEta_->Fill(eta, error);
0381 IPErrVsPt_->Fill(pt, error);
0382 IPErrVsEtaVsPhi_->Fill(eta, phi, error);
0383 }
0384
0385
0386
0387
0388 GeneralPurposeVertexAnalyzer::GeneralPurposeVertexAnalyzer(const edm::ParameterSet &iConfig)
0389 : conf_(iConfig),
0390 ndof_(iConfig.getParameter<int>("ndof")),
0391 errorPrinted_(false),
0392 vertexInputTag_(iConfig.getParameter<edm::InputTag>("vertexLabel")),
0393 beamSpotInputTag_(iConfig.getParameter<edm::InputTag>("beamSpotLabel")),
0394 vertexToken_(consumes<reco::VertexCollection>(vertexInputTag_)),
0395 scoreToken_(consumes<VertexScore>(vertexInputTag_)),
0396 beamspotToken_(consumes<reco::BeamSpot>(beamSpotInputTag_)),
0397
0398 vposx_(iConfig.getParameter<double>("Xpos")),
0399 vposy_(iConfig.getParameter<double>("Ypos")),
0400 tkNoBin_(iConfig.getParameter<int>("TkSizeBin")),
0401 tkNoMin_(iConfig.getParameter<double>("TkSizeMin")),
0402 tkNoMax_(iConfig.getParameter<double>("TkSizeMax")),
0403 dxyBin_(iConfig.getParameter<int>("DxyBin")),
0404 dxyMin_(iConfig.getParameter<double>("DxyMin")),
0405 dxyMax_(iConfig.getParameter<double>("DxyMax")),
0406 dzBin_(iConfig.getParameter<int>("DzBin")),
0407 dzMin_(iConfig.getParameter<double>("DzMin")),
0408 dzMax_(iConfig.getParameter<double>("DzMax")),
0409 phiBin_(iConfig.getParameter<int>("PhiBin")),
0410 phiBin2D_(iConfig.getParameter<int>("PhiBin2D")),
0411 phiMin_(iConfig.getParameter<double>("PhiMin")),
0412 phiMax_(iConfig.getParameter<double>("PhiMax")),
0413 etaBin_(iConfig.getParameter<int>("EtaBin")),
0414 etaBin2D_(iConfig.getParameter<int>("EtaBin2D")),
0415 etaMin_(iConfig.getParameter<double>("EtaMin")),
0416 etaMax_(iConfig.getParameter<double>("EtaMax")),
0417
0418 nbvtx(nullptr),
0419 bsX(nullptr),
0420 bsY(nullptr),
0421 bsZ(nullptr),
0422 bsSigmaZ(nullptr),
0423 bsDxdz(nullptr),
0424 bsDydz(nullptr),
0425 bsBeamWidthX(nullptr),
0426 bsBeamWidthY(nullptr),
0427 bsType(nullptr),
0428 trackpt(nullptr),
0429 sumpt(nullptr),
0430 ntracks(nullptr),
0431 weight(nullptr),
0432 chi2ndf(nullptr),
0433 chi2prob(nullptr),
0434 dxy2(nullptr),
0435 phi_pt1(nullptr),
0436 eta_pt1(nullptr),
0437 phi_pt10(nullptr),
0438 eta_pt10(nullptr) {
0439 usesResource(TFileService::kSharedResource);
0440 }
0441
0442
0443
0444
0445
0446
0447 void GeneralPurposeVertexAnalyzer::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0448 using namespace edm;
0449
0450 const auto &recVtxs = iEvent.getHandle(vertexToken_);
0451 const auto &scores = iEvent.getHandle(scoreToken_);
0452 const auto &beamSpotHandle = iEvent.getHandle(beamspotToken_);
0453 reco::BeamSpot beamSpot = *beamSpotHandle;
0454
0455
0456 if (recVtxs.isValid() == false || beamSpotHandle.isValid() == false) {
0457 edm::LogWarning("GeneralPurposeVertexAnalyzer")
0458 << " Some products not available in the event: VertexCollection " << vertexInputTag_ << " " << recVtxs.isValid()
0459 << " BeamSpot " << beamSpotInputTag_ << " " << beamSpotHandle.isValid() << ". Skipping plots for this event";
0460 return;
0461 }
0462
0463
0464 bool ok{true};
0465 for (const auto &v : *recVtxs) {
0466 if (v.tracksSize() > 0) {
0467 const auto &ref = v.trackRefAt(0);
0468 if (ref.isNull() || !ref.isAvailable()) {
0469 if (!errorPrinted_) {
0470 edm::LogWarning("GeneralPurposeVertexAnalyzer")
0471 << "Skipping vertex collection: " << vertexInputTag_
0472 << " since likely the track collection the vertex has refs pointing to is missing (at least the first "
0473 "TrackBaseRef is null or not available)";
0474 } else {
0475 errorPrinted_ = true;
0476 }
0477 ok = false;
0478 }
0479 }
0480 }
0481
0482 if (!ok) {
0483 return;
0484 }
0485
0486 nbvtx->Fill(recVtxs->size());
0487 int ng = 0;
0488 for (auto const &vx : (*recVtxs)) {
0489 if (vx.isValid() && !vx.isFake() && vx.ndof() >= ndof_) {
0490 ++ng;
0491 }
0492 }
0493 nbgvtx->Fill(ng);
0494
0495 if (scores.isValid() && !(*scores).empty()) {
0496 auto pvScore = (*scores).get(0);
0497 score[1]->Fill(std::sqrt(pvScore));
0498 for (unsigned int i = 1; i < (*scores).size(); ++i) {
0499 score[0]->Fill(std::sqrt((*scores).get(i)));
0500 }
0501 }
0502
0503
0504 if (!recVtxs->empty()) {
0505 vertexPlots(recVtxs->front(), beamSpot, 1);
0506 pvTracksPlots(recVtxs->front());
0507
0508 for (reco::VertexCollection::const_iterator v = recVtxs->begin() + 1; v != recVtxs->end(); ++v) {
0509 vertexPlots(*v, beamSpot, 0);
0510 }
0511 }
0512
0513
0514 bsX->Fill(beamSpot.x0());
0515 bsY->Fill(beamSpot.y0());
0516 bsZ->Fill(beamSpot.z0());
0517 bsSigmaZ->Fill(beamSpot.sigmaZ());
0518 bsDxdz->Fill(beamSpot.dxdz());
0519 bsDydz->Fill(beamSpot.dydz());
0520 bsBeamWidthX->Fill(beamSpot.BeamWidthX() * cmToUm);
0521 bsBeamWidthY->Fill(beamSpot.BeamWidthY() * cmToUm);
0522 bsType->Fill(beamSpot.type());
0523 }
0524
0525 void GeneralPurposeVertexAnalyzer::pvTracksPlots(const reco::Vertex &v) {
0526 if (!v.isValid())
0527 return;
0528 if (v.isFake())
0529 return;
0530
0531 if (v.tracksSize() == 0) {
0532 ntracks->Fill(0);
0533 return;
0534 }
0535
0536 const math::XYZPoint myVertex(v.position().x(), v.position().y(), v.position().z());
0537
0538 float sumPT = 0.;
0539 for (const auto &t : v.tracks()) {
0540 const bool isHighPurity = t->quality(reco::TrackBase::highPurity);
0541 if (!isHighPurity) {
0542 continue;
0543 }
0544
0545 const float pt = t->pt();
0546 trackpt->Fill(pt);
0547
0548 if (pt < 1.f) {
0549 continue;
0550 }
0551
0552 const float pt2 = pt * pt;
0553 const float eta = t->eta();
0554 const float phi = t->phi();
0555 const float w = v.trackWeight(t);
0556 const float chi2NDF = t->normalizedChi2();
0557 const float chi2Prob = TMath::Prob(t->chi2(), static_cast<int>(t->ndof()));
0558 const float Dxy = t->dxy(myVertex) * cmToUm;
0559
0560 sumPT += pt2;
0561
0562 weight->Fill(w);
0563 chi2ndf->Fill(chi2NDF);
0564 chi2prob->Fill(chi2Prob);
0565 dxy2->Fill(Dxy);
0566 phi_pt1->Fill(phi);
0567 eta_pt1->Fill(eta);
0568
0569
0570 dxy_pt1.fillIPMonitor(t, myVertex);
0571
0572
0573 dz_pt1.fillIPMonitor(t, myVertex);
0574
0575 if (pt >= 10.f) {
0576 phi_pt10->Fill(phi);
0577 eta_pt10->Fill(eta);
0578
0579
0580 dxy_pt10.fillIPMonitor(t, myVertex);
0581
0582
0583 dz_pt10.fillIPMonitor(t, myVertex);
0584 }
0585 }
0586
0587 ntracks->Fill(static_cast<float>(v.tracks().size()));
0588 sumpt->Fill(sumPT);
0589 }
0590
0591 void GeneralPurposeVertexAnalyzer::vertexPlots(const reco::Vertex &v, const reco::BeamSpot &beamSpot, int i) {
0592 if (i < 0 || i > 1)
0593 return;
0594 if (!v.isValid())
0595 type[i]->Fill(2.);
0596 else if (v.isFake())
0597 type[i]->Fill(1.);
0598 else
0599 type[i]->Fill(0.);
0600
0601 if (v.isValid() && !v.isFake()) {
0602 float weight = 0;
0603 for (reco::Vertex::trackRef_iterator t = v.tracks_begin(); t != v.tracks_end(); t++)
0604 weight += v.trackWeight(*t);
0605 trksWeight[i]->Fill(weight);
0606 nbtksinvtx[i]->Fill(v.tracksSize());
0607 ntracksVsZ[i]->Fill(v.position().z() - beamSpot.z0(), v.tracksSize());
0608
0609 vtxchi2[i]->Fill(v.chi2());
0610 vtxndf[i]->Fill(v.ndof());
0611 vtxprob[i]->Fill(ChiSquaredProbability(v.chi2(), v.ndof()));
0612
0613 xrec[i]->Fill(v.position().x());
0614 yrec[i]->Fill(v.position().y());
0615 zrec[i]->Fill(v.position().z());
0616
0617 float xb = beamSpot.x0() + beamSpot.dxdz() * (v.position().z() - beamSpot.z0());
0618 float yb = beamSpot.y0() + beamSpot.dydz() * (v.position().z() - beamSpot.z0());
0619 xDiff[i]->Fill((v.position().x() - xb) * cmToUm);
0620 yDiff[i]->Fill((v.position().y() - yb) * cmToUm);
0621
0622 xerr[i]->Fill(v.xError() * cmToUm);
0623 yerr[i]->Fill(v.yError() * cmToUm);
0624 zerr[i]->Fill(v.zError() * cmToUm);
0625 xerrVsTrks[i]->Fill(weight, v.xError() * cmToUm);
0626 yerrVsTrks[i]->Fill(weight, v.yError() * cmToUm);
0627 zerrVsTrks[i]->Fill(weight, v.zError() * cmToUm);
0628
0629 nans[i]->Fill(1., edm::isNotFinite(v.position().x()) * 1.);
0630 nans[i]->Fill(2., edm::isNotFinite(v.position().y()) * 1.);
0631 nans[i]->Fill(3., edm::isNotFinite(v.position().z()) * 1.);
0632
0633 int index = 3;
0634 for (int k = 0; k != 3; k++) {
0635 for (int j = k; j != 3; j++) {
0636 index++;
0637 nans[i]->Fill(index * 1., edm::isNotFinite(v.covariance(k, j)) * 1.);
0638
0639 if (j == k && v.covariance(k, j) < 0) {
0640 nans[i]->Fill(index * 1., 1.);
0641 }
0642 }
0643 }
0644 }
0645 }
0646
0647 template <typename T, typename... Args>
0648 T *GeneralPurposeVertexAnalyzer::book(const Args &...args) const {
0649 T *t = fs_->make<T>(args...);
0650 return t;
0651 }
0652
0653
0654 void GeneralPurposeVertexAnalyzer::beginJob() {
0655 nbvtx = book<TH1I>("vtxNbr", "Reconstructed Vertices in Event", 80, -0.5, 79.5);
0656 nbgvtx = book<TH1I>("goodvtxNbr", "Reconstructed Good Vertices in Event", 80, -0.5, 79.5);
0657
0658 nbtksinvtx[0] = book<TH1D>("otherVtxTrksNbr", "Reconstructed Tracks in Vertex (other Vtx)", 40, -0.5, 99.5);
0659 ntracksVsZ[0] =
0660 book<TProfile>("otherVtxTrksVsZ", "Reconstructed Tracks in Vertex (other Vtx) vs Z", 80, -20., 20., 0., 100., "");
0661 ntracksVsZ[0]->SetXTitle("z-bs");
0662 ntracksVsZ[0]->SetYTitle("#tracks");
0663
0664 score[0] = book<TH1D>("otherVtxScore", "sqrt(score) (other Vtx)", 100, 0., 400.);
0665 trksWeight[0] = book<TH1D>("otherVtxTrksWeight", "Total weight of Tracks in Vertex (other Vtx)", 40, 0, 100.);
0666 vtxchi2[0] = book<TH1D>("otherVtxChi2", "#chi^{2} (other Vtx)", 100, 0., 200.);
0667 vtxndf[0] = book<TH1D>("otherVtxNdf", "ndof (other Vtx)", 100, 0., 200.);
0668 vtxprob[0] = book<TH1D>("otherVtxProb", "#chi^{2} probability (other Vtx)", 100, 0., 1.);
0669 nans[0] = book<TH1D>("otherVtxNans", "Illegal values for x,y,z,xx,xy,xz,yy,yz,zz (other Vtx)", 9, 0.5, 9.5);
0670
0671 nbtksinvtx[1] = book<TH1D>("tagVtxTrksNbr", "Reconstructed Tracks in Vertex (tagged Vtx)", 100, -0.5, 99.5);
0672 ntracksVsZ[1] =
0673 book<TProfile>("tagVtxTrksVsZ", "Reconstructed Tracks in Vertex (tagged Vtx) vs Z", 80, -20., 20., 0., 100., "");
0674 ntracksVsZ[1]->SetXTitle("z-bs");
0675 ntracksVsZ[1]->SetYTitle("#tracks");
0676
0677 score[1] = book<TH1D>("tagVtxScore", "sqrt(score) (tagged Vtx)", 100, 0., 400.);
0678 trksWeight[1] = book<TH1D>("tagVtxTrksWeight", "Total weight of Tracks in Vertex (tagged Vtx)", 100, 0, 100.);
0679 vtxchi2[1] = book<TH1D>("tagVtxChi2", "#chi^{2} (tagged Vtx)", 100, 0., 200.);
0680 vtxndf[1] = book<TH1D>("tagVtxNdf", "ndof (tagged Vtx)", 100, 0., 200.);
0681 vtxprob[1] = book<TH1D>("tagVtxProb", "#chi^{2} probability (tagged Vtx)", 100, 0., 1.);
0682 nans[1] = book<TH1D>("tagVtxNans", "Illegal values for x,y,z,xx,xy,xz,yy,yz,zz (tagged Vtx)", 9, 0.5, 9.5);
0683
0684 xrec[0] = book<TH1D>("otherPosX", "Position x Coordinate (other Vtx)", 100, vposx_ - 0.1, vposx_ + 0.1);
0685 yrec[0] = book<TH1D>("otherPosY", "Position y Coordinate (other Vtx)", 100, vposy_ - 0.1, vposy_ + 0.1);
0686 zrec[0] = book<TH1D>("otherPosZ", "Position z Coordinate (other Vtx)", 100, -20., 20.);
0687 xDiff[0] = book<TH1D>("otherDiffX", "X distance from BeamSpot (other Vtx)", 100, -500, 500);
0688 yDiff[0] = book<TH1D>("otherDiffY", "Y distance from BeamSpot (other Vtx)", 100, -500, 500);
0689 xerr[0] = book<TH1D>("otherErrX", "Uncertainty x Coordinate (other Vtx)", 100, 0., 100);
0690 yerr[0] = book<TH1D>("otherErrY", "Uncertainty y Coordinate (other Vtx)", 100, 0., 100);
0691 zerr[0] = book<TH1D>("otherErrZ", "Uncertainty z Coordinate (other Vtx)", 100, 0., 100);
0692 xerrVsTrks[0] = book<TH2D>(
0693 "otherErrVsWeightX", "Uncertainty x Coordinate vs. track weight (other Vtx)", 100, 0, 100., 100, 0., 100);
0694 yerrVsTrks[0] = book<TH2D>(
0695 "otherErrVsWeightY", "Uncertainty y Coordinate vs. track weight (other Vtx)", 100, 0, 100., 100, 0., 100);
0696 zerrVsTrks[0] = book<TH2D>(
0697 "otherErrVsWeightZ", "Uncertainty z Coordinate vs. track weight (other Vtx)", 100, 0, 100., 100, 0., 100);
0698
0699 xrec[1] = book<TH1D>("tagPosX", "Position x Coordinate (tagged Vtx)", 100, vposx_ - 0.1, vposx_ + 0.1);
0700 yrec[1] = book<TH1D>("tagPosY", "Position y Coordinate (tagged Vtx)", 100, vposy_ - 0.1, vposy_ + 0.1);
0701 zrec[1] = book<TH1D>("tagPosZ", "Position z Coordinate (tagged Vtx)", 100, -20., 20.);
0702 xDiff[1] = book<TH1D>("tagDiffX", "X distance from BeamSpot (tagged Vtx)", 100, -500, 500);
0703 yDiff[1] = book<TH1D>("tagDiffY", "Y distance from BeamSpot (tagged Vtx)", 100, -500, 500);
0704 xerr[1] = book<TH1D>("tagErrX", "Uncertainty x Coordinate (tagged Vtx)", 100, 0., 100);
0705 yerr[1] = book<TH1D>("tagErrY", "Uncertainty y Coordinate (tagged Vtx)", 100, 0., 100);
0706 zerr[1] = book<TH1D>("tagErrZ", "Uncertainty z Coordinate (tagged Vtx)", 100, 0., 100);
0707 xerrVsTrks[1] = book<TH2D>(
0708 "tagErrVsWeightX", "Uncertainty x Coordinate vs. track weight (tagged Vtx)", 100, 0, 100., 100, 0., 100);
0709 yerrVsTrks[1] = book<TH2D>(
0710 "tagErrVsWeightY", "Uncertainty y Coordinate vs. track weight (tagged Vtx)", 100, 0, 100., 100, 0., 100);
0711 zerrVsTrks[1] = book<TH2D>(
0712 "tagErrVsWeightZ", "Uncertainty z Coordinate vs. track weight (tagged Vtx)", 100, 0, 100., 100, 0., 100);
0713
0714 type[0] = book<TH1D>("otherType", "Vertex type (other Vtx)", 3, -0.5, 2.5);
0715 type[1] = book<TH1D>("tagType", "Vertex type (tagged Vtx)", 3, -0.5, 2.5);
0716
0717 for (int i = 0; i < 2; ++i) {
0718 type[i]->GetXaxis()->SetBinLabel(1, "Valid, real");
0719 type[i]->GetXaxis()->SetBinLabel(2, "Valid, fake");
0720 type[i]->GetXaxis()->SetBinLabel(3, "Invalid");
0721 }
0722
0723 bsX = book<TH1D>("bsX", "BeamSpot x0", 100, vposx_ - 0.1, vposx_ + 0.1);
0724 bsY = book<TH1D>("bsY", "BeamSpot y0", 100, vposy_ - 0.1, vposy_ + 0.1);
0725 bsZ = book<TH1D>("bsZ", "BeamSpot z0", 100, -2., 2.);
0726 bsSigmaZ = book<TH1D>("bsSigmaZ", "BeamSpot sigmaZ", 100, 0., 10.);
0727 bsDxdz = book<TH1D>("bsDxdz", "BeamSpot dxdz", 100, -0.0003, 0.0003);
0728 bsDydz = book<TH1D>("bsDydz", "BeamSpot dydz", 100, -0.0003, 0.0003);
0729 bsBeamWidthX = book<TH1D>("bsBeamWidthX", "BeamSpot BeamWidthX", 100, 0., 100.);
0730 bsBeamWidthY = book<TH1D>("bsBeamWidthY", "BeamSpot BeamWidthY", 100, 0., 100.);
0731 bsType = book<TH1D>("bsType", "BeamSpot type", 4, -1.5, 2.5);
0732 bsType->GetXaxis()->SetBinLabel(1, "Unknown");
0733 bsType->GetXaxis()->SetBinLabel(2, "Fake");
0734 bsType->GetXaxis()->SetBinLabel(3, "LHC");
0735 bsType->GetXaxis()->SetBinLabel(4, "Tracker");
0736
0737
0738 std::string s_1 = "PV Tracks (p_{T} > 1 GeV)";
0739 std::string s_10 = "PV Tracks (p_{T} > 10 GeV)";
0740
0741 ntracks = book<TH1D>("ntracks", fmt::sprintf("number of %s", s_1).c_str(), tkNoBin_, tkNoMin_, tkNoMax_);
0742 ntracks->SetXTitle(fmt::sprintf("Number of %s per Event", s_1).c_str());
0743 ntracks->SetYTitle("Number of Events");
0744
0745 weight = book<TH1D>("weight", fmt::sprintf("weight of %s", s_1).c_str(), 100, 0., 1.);
0746 weight->SetXTitle(fmt::sprintf("weight of %s per Event", s_1).c_str());
0747 weight->SetYTitle("Number of Events");
0748
0749 sumpt = book<TH1D>("sumpt", fmt::sprintf("#Sum p_{T} of %s", s_1).c_str(), 100, -0.5, 249.5);
0750 chi2ndf = book<TH1D>("chi2ndf", fmt::sprintf("%s #chi^{2}/ndof", s_1).c_str(), 100, 0., 20.);
0751 chi2prob = book<TH1D>("chi2prob", fmt::sprintf("%s #chi^{2} probability", s_1).c_str(), 100, 0., 1.);
0752
0753 dxy2 = book<TH1D>("dxyzoom", fmt::sprintf("%s d_{xy} (#mum)", s_1).c_str(), dxyBin_, dxyMin_ / 5., dxyMax_ / 5.);
0754
0755 trackpt = gpVertexAnalyzer::makeTH1IfLog(
0756 fs_, true, false, "pt_track", "PV tracks p_{T};PV tracks p_{T} [GeV];#tracks", 49, log10(1.), log10(50.));
0757 phi_pt1 =
0758 book<TH1D>("phi_pt1", fmt::sprintf("%s #phi; PV tracks #phi;#tracks", s_1).c_str(), phiBin_, phiMin_, phiMax_);
0759 eta_pt1 =
0760 book<TH1D>("eta_pt1", fmt::sprintf("%s #eta; PV tracks #eta;#tracks", s_1).c_str(), etaBin_, etaMin_, etaMax_);
0761 phi_pt10 =
0762 book<TH1D>("phi_pt10", fmt::sprintf("%s #phi; PV tracks #phi;#tracks", s_10).c_str(), phiBin_, phiMin_, phiMax_);
0763 eta_pt10 =
0764 book<TH1D>("eta_pt10", fmt::sprintf("%s #phi; PV tracks #eta;#tracks", s_10).c_str(), etaBin_, etaMin_, etaMax_);
0765
0766
0767 dxy_pt1.varname_ = "xy";
0768 dxy_pt1.pTcut_ = 1.f;
0769 dxy_pt1.bookIPMonitor(conf_, fs_);
0770
0771 dxy_pt10.varname_ = "xy";
0772 dxy_pt10.pTcut_ = 10.f;
0773 dxy_pt10.bookIPMonitor(conf_, fs_);
0774
0775 dz_pt1.varname_ = "z";
0776 dz_pt1.pTcut_ = 1.f;
0777 dz_pt1.bookIPMonitor(conf_, fs_);
0778
0779 dz_pt10.varname_ = "z";
0780 dz_pt10.pTcut_ = 10.f;
0781 dz_pt10.bookIPMonitor(conf_, fs_);
0782 }
0783
0784
0785 void GeneralPurposeVertexAnalyzer::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0786 edm::ParameterSetDescription desc;
0787 desc.add<int>("ndof", 4);
0788 desc.add<edm::InputTag>("vertexLabel", edm::InputTag("offlinePrimaryVertices"));
0789 desc.add<edm::InputTag>("beamSpotLabel", edm::InputTag("offlineBeamSpot"));
0790 desc.add<double>("Xpos", 0.1);
0791 desc.add<double>("Ypos", 0.0);
0792 desc.add<int>("TkSizeBin", 100);
0793 desc.add<double>("TkSizeMin", 499.5);
0794 desc.add<double>("TkSizeMax", -0.5);
0795 desc.add<int>("DxyBin", 100);
0796 desc.add<double>("DxyMin", -2000.);
0797 desc.add<double>("DxyMax", 2000.);
0798 desc.add<int>("DzBin", 100);
0799 desc.add<double>("DzMin", -2000.0);
0800 desc.add<double>("DzMax", 2000.0);
0801 desc.add<int>("PhiBin", 32);
0802 desc.add<int>("PhiBin2D", 12);
0803 desc.add<double>("PhiMin", -M_PI);
0804 desc.add<double>("PhiMax", M_PI);
0805 desc.add<int>("EtaBin", 26);
0806 desc.add<int>("EtaBin2D", 8);
0807 desc.add<double>("EtaMin", -2.7);
0808 desc.add<double>("EtaMax", 2.7);
0809 desc.add<int>("PtBin", 49);
0810 desc.add<double>("PtMin", 1.);
0811 desc.add<double>("PtMax", 50.);
0812 descriptions.addWithDefaultLabel(desc);
0813 }
0814
0815
0816 DEFINE_FWK_MODULE(GeneralPurposeVertexAnalyzer);
