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