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