Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-11-17 23:29:39

0001 // user includes
0002 #include "CommonTools/Statistics/interface/ChiSquaredProbability.h"
0003 #include "DQMOffline/RecoB/plugins/PrimaryVertexMonitor.h"
0004 #include "FWCore/Framework/interface/MakerMacros.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 #include "FWCore/ServiceRegistry/interface/Service.h"
0007 #include "FWCore/Utilities/interface/isFinite.h"
0008 
0009 // ROOT includes
0010 #include "TMath.h"
0011 
0012 // system includes
0013 #include <fmt/format.h>
0014 
0015 using namespace reco;
0016 using namespace edm;
0017 
0018 PrimaryVertexMonitor::PrimaryVertexMonitor(const edm::ParameterSet& pSet)
0019     : vertexInputTag_(pSet.getParameter<InputTag>("vertexLabel")),
0020       beamSpotInputTag_(pSet.getParameter<InputTag>("beamSpotLabel")),
0021       vertexToken_(consumes<reco::VertexCollection>(vertexInputTag_)),
0022       scoreToken_(consumes<VertexScore>(vertexInputTag_)),
0023       beamspotToken_(consumes<reco::BeamSpot>(beamSpotInputTag_)),
0024       conf_(pSet),
0025       TopFolderName_(pSet.getParameter<std::string>("TopFolderName")),
0026       AlignmentLabel_(pSet.getParameter<std::string>("AlignmentLabel")),
0027       ndof_(pSet.getParameter<int>("ndof")),
0028       useHPfoAlignmentPlots_(pSet.getParameter<bool>("useHPforAlignmentPlots")),
0029       errorPrinted_(false),
0030       nbvtx(nullptr),
0031       bsX(nullptr),
0032       bsY(nullptr),
0033       bsZ(nullptr),
0034       bsSigmaZ(nullptr),
0035       bsDxdz(nullptr),
0036       bsDydz(nullptr),
0037       bsBeamWidthX(nullptr),
0038       bsBeamWidthY(nullptr),
0039       bsType(nullptr),
0040       sumpt(nullptr),
0041       ntracks(nullptr),
0042       weight(nullptr),
0043       chi2ndf(nullptr),
0044       chi2prob(nullptr),
0045       trackpt(nullptr),
0046       phi_pt1(nullptr),
0047       eta_pt1(nullptr),
0048       phi_pt10(nullptr),
0049       eta_pt10(nullptr),
0050       dxy2(nullptr) {}
0051 
0052 // -- BeginRun
0053 //---------------------------------------------------------------------------------//
0054 void PrimaryVertexMonitor::bookHistograms(DQMStore::IBooker& iBooker, edm::Run const&, edm::EventSetup const&) {
0055   std::string dqmLabel = "";
0056 
0057   //
0058   // Book all histograms.
0059   //
0060 
0061   //  get the store
0062   dqmLabel = TopFolderName_ + "/" + vertexInputTag_.label();
0063   iBooker.setCurrentFolder(dqmLabel);
0064 
0065   //   xPos = iBooker.book1D ("xPos","x Coordinate" ,100, -0.1, 0.1);
0066   auto maxPU = conf_.getParameter<double>("PUMax");
0067   nbvtx = iBooker.book1D("vtxNbr", "Reconstructed Vertices in Event", maxPU, -0.5, maxPU - 0.5);
0068   nbgvtx = iBooker.book1D("goodvtxNbr", "Reconstructed Good Vertices in Event", maxPU, -0.5, maxPU - 0.5);
0069 
0070   // to be configured each year...
0071   auto vposx = conf_.getParameter<double>("Xpos");
0072   auto vposy = conf_.getParameter<double>("Ypos");
0073 
0074   nbtksinvtx[0] = iBooker.book1D("otherVtxTrksNbr", "Reconstructed Tracks in Vertex (other Vtx)", 40, -0.5, 99.5);
0075   ntracksVsZ[0] = iBooker.bookProfile(
0076       "otherVtxTrksVsZ", "Reconstructed Tracks in Vertex (other Vtx) vs Z", 80, -20., 20., 50, 0, 100, "");
0077   ntracksVsZ[0]->setAxisTitle("z-bs", 1);
0078   ntracksVsZ[0]->setAxisTitle("#tracks", 2);
0079 
0080   score[0] = iBooker.book1D("otherVtxScore", "sqrt(score) (other Vtx)", 100, 0., 400.);
0081   trksWeight[0] = iBooker.book1D("otherVtxTrksWeight", "Total weight of Tracks in Vertex (other Vtx)", 40, 0, 100.);
0082   vtxchi2[0] = iBooker.book1D("otherVtxChi2", "#chi^{2} (other Vtx)", 100, 0., 200.);
0083   vtxndf[0] = iBooker.book1D("otherVtxNdf", "ndof (other Vtx)", 100, 0., 200.);
0084   vtxprob[0] = iBooker.book1D("otherVtxProb", "#chi^{2} probability (other Vtx)", 100, 0., 1.);
0085   nans[0] = iBooker.book1D("otherVtxNans", "Illegal values for x,y,z,xx,xy,xz,yy,yz,zz (other Vtx)", 9, 0.5, 9.5);
0086 
0087   nbtksinvtx[1] = iBooker.book1D("tagVtxTrksNbr", "Reconstructed Tracks in Vertex (tagged Vtx)", 100, -0.5, 99.5);
0088   ntracksVsZ[1] = iBooker.bookProfile(
0089       "tagVtxTrksVsZ", "Reconstructed Tracks in Vertex (tagged Vtx) vs Z", 80, -20., 20., 50, 0, 100, "");
0090   ntracksVsZ[1]->setAxisTitle("z-bs", 1);
0091   ntracksVsZ[1]->setAxisTitle("#tracks", 2);
0092 
0093   score[1] = iBooker.book1D("tagVtxScore", "sqrt(score) (tagged Vtx)", 100, 0., 400.);
0094   trksWeight[1] = iBooker.book1D("tagVtxTrksWeight", "Total weight of Tracks in Vertex (tagged Vtx)", 100, 0, 100.);
0095   vtxchi2[1] = iBooker.book1D("tagVtxChi2", "#chi^{2} (tagged Vtx)", 100, 0., 200.);
0096   vtxndf[1] = iBooker.book1D("tagVtxNdf", "ndof (tagged Vtx)", 100, 0., 200.);
0097   vtxprob[1] = iBooker.book1D("tagVtxProb", "#chi^{2} probability (tagged Vtx)", 100, 0., 1.);
0098   nans[1] = iBooker.book1D("tagVtxNans", "Illegal values for x,y,z,xx,xy,xz,yy,yz,zz (tagged Vtx)", 9, 0.5, 9.5);
0099 
0100   xrec[0] = iBooker.book1D("otherPosX", "Position x Coordinate (other Vtx)", 100, vposx - 0.1, vposx + 0.1);
0101   yrec[0] = iBooker.book1D("otherPosY", "Position y Coordinate (other Vtx)", 100, vposy - 0.1, vposy + 0.1);
0102   zrec[0] = iBooker.book1D("otherPosZ", "Position z Coordinate (other Vtx)", 100, -20., 20.);
0103   xDiff[0] = iBooker.book1D("otherDiffX", "X distance from BeamSpot (other Vtx)", 100, -500, 500);
0104   yDiff[0] = iBooker.book1D("otherDiffY", "Y distance from BeamSpot (other Vtx)", 100, -500, 500);
0105   xerr[0] = iBooker.book1D("otherErrX", "Uncertainty x Coordinate (other Vtx)", 100, 0., 100);
0106   yerr[0] = iBooker.book1D("otherErrY", "Uncertainty y Coordinate (other Vtx)", 100, 0., 100);
0107   zerr[0] = iBooker.book1D("otherErrZ", "Uncertainty z Coordinate (other Vtx)", 100, 0., 100);
0108   xerrVsTrks[0] = iBooker.book2D(
0109       "otherErrVsWeightX", "Uncertainty x Coordinate vs. track weight (other Vtx)", 100, 0, 100., 100, 0., 100);
0110   yerrVsTrks[0] = iBooker.book2D(
0111       "otherErrVsWeightY", "Uncertainty y Coordinate vs. track weight (other Vtx)", 100, 0, 100., 100, 0., 100);
0112   zerrVsTrks[0] = iBooker.book2D(
0113       "otherErrVsWeightZ", "Uncertainty z Coordinate vs. track weight (other Vtx)", 100, 0, 100., 100, 0., 100);
0114 
0115   xrec[1] = iBooker.book1D("tagPosX", "Position x Coordinate (tagged Vtx)", 100, vposx - 0.1, vposx + 0.1);
0116   yrec[1] = iBooker.book1D("tagPosY", "Position y Coordinate (tagged Vtx)", 100, vposy - 0.1, vposy + 0.1);
0117   zrec[1] = iBooker.book1D("tagPosZ", "Position z Coordinate (tagged Vtx)", 100, -20., 20.);
0118   xDiff[1] = iBooker.book1D("tagDiffX", "X distance from BeamSpot (tagged Vtx)", 100, -500, 500);
0119   yDiff[1] = iBooker.book1D("tagDiffY", "Y distance from BeamSpot (tagged Vtx)", 100, -500, 500);
0120   xerr[1] = iBooker.book1D("tagErrX", "Uncertainty x Coordinate (tagged Vtx)", 100, 0., 100);
0121   yerr[1] = iBooker.book1D("tagErrY", "Uncertainty y Coordinate (tagged Vtx)", 100, 0., 100);
0122   zerr[1] = iBooker.book1D("tagErrZ", "Uncertainty z Coordinate (tagged Vtx)", 100, 0., 100);
0123   xerrVsTrks[1] = iBooker.book2D(
0124       "tagErrVsWeightX", "Uncertainty x Coordinate vs. track weight (tagged Vtx)", 100, 0, 100., 100, 0., 100);
0125   yerrVsTrks[1] = iBooker.book2D(
0126       "tagErrVsWeightY", "Uncertainty y Coordinate vs. track weight (tagged Vtx)", 100, 0, 100., 100, 0., 100);
0127   zerrVsTrks[1] = iBooker.book2D(
0128       "tagErrVsWeightZ", "Uncertainty z Coordinate vs. track weight (tagged Vtx)", 100, 0, 100., 100, 0., 100);
0129 
0130   type[0] = iBooker.book1D("otherType", "Vertex type (other Vtx)", 3, -0.5, 2.5);
0131   type[1] = iBooker.book1D("tagType", "Vertex type (tagged Vtx)", 3, -0.5, 2.5);
0132   for (int i = 0; i < 2; ++i) {
0133     type[i]->setBinLabel(1, "Valid, real");
0134     type[i]->setBinLabel(2, "Valid, fake");
0135     type[i]->setBinLabel(3, "Invalid");
0136   }
0137 
0138   //  get the store
0139   dqmLabel = TopFolderName_ + "/" + beamSpotInputTag_.label();
0140   iBooker.setCurrentFolder(dqmLabel);
0141 
0142   bsX = iBooker.book1D("bsX", "BeamSpot x0", 100, vposx - 0.1, vposx + 0.1);
0143   bsY = iBooker.book1D("bsY", "BeamSpot y0", 100, vposy - 0.1, vposy + 0.1);
0144   bsZ = iBooker.book1D("bsZ", "BeamSpot z0", 100, -2., 2.);
0145   bsSigmaZ = iBooker.book1D("bsSigmaZ", "BeamSpot sigmaZ", 100, 0., 10.);
0146   bsDxdz = iBooker.book1D("bsDxdz", "BeamSpot dxdz", 100, -0.0003, 0.0003);
0147   bsDydz = iBooker.book1D("bsDydz", "BeamSpot dydz", 100, -0.0003, 0.0003);
0148   bsBeamWidthX = iBooker.book1D("bsBeamWidthX", "BeamSpot BeamWidthX", 500, 0., 15.);
0149   bsBeamWidthY = iBooker.book1D("bsBeamWidthY", "BeamSpot BeamWidthY", 500, 0., 15.);
0150   bsType = iBooker.book1D("bsType", "BeamSpot type", 4, -1.5, 2.5);
0151   bsType->setBinLabel(1, "Unknown");
0152   bsType->setBinLabel(2, "Fake");
0153   bsType->setBinLabel(3, "LHC");
0154   bsType->setBinLabel(4, "Tracker");
0155 
0156   //  get the store
0157   dqmLabel = TopFolderName_ + "/" + AlignmentLabel_;
0158   iBooker.setCurrentFolder(dqmLabel);
0159 
0160   int TKNoBin = conf_.getParameter<int>("TkSizeBin");
0161   double TKNoMin = conf_.getParameter<double>("TkSizeMin");
0162   double TKNoMax = conf_.getParameter<double>("TkSizeMax");
0163 
0164   int DxyBin = conf_.getParameter<int>("DxyBin");
0165   double DxyMin = conf_.getParameter<double>("DxyMin");
0166   double DxyMax = conf_.getParameter<double>("DxyMax");
0167 
0168   int PhiBin = conf_.getParameter<int>("PhiBin");
0169   double PhiMin = conf_.getParameter<double>("PhiMin");
0170   double PhiMax = conf_.getParameter<double>("PhiMax");
0171 
0172   int EtaBin = conf_.getParameter<int>("EtaBin");
0173   double EtaMin = conf_.getParameter<double>("EtaMin");
0174   double EtaMax = conf_.getParameter<double>("EtaMax");
0175 
0176   int PtBin = conf_.getParameter<int>("PtBin");
0177   double PtMin = conf_.getParameter<double>("PtMin");
0178   double PtMax = conf_.getParameter<double>("PtMax");
0179 
0180   ntracks = iBooker.book1D("ntracks", "number of PV tracks (p_{T} > 1 GeV)", TKNoBin, TKNoMin, TKNoMax);
0181   ntracks->setAxisTitle("Number of PV Tracks (p_{T} > 1 GeV) per Event", 1);
0182   ntracks->setAxisTitle("Number of Event", 2);
0183 
0184   weight = iBooker.book1D("weight", "weight of PV tracks (p_{T} > 1 GeV)", 100, 0., 1.);
0185   weight->setAxisTitle("weight of PV Tracks (p_{T} > 1 GeV) per Event", 1);
0186   weight->setAxisTitle("Number of Event", 2);
0187 
0188   sumpt = iBooker.book1D("sumpt", "#Sum p_{T} of PV tracks (p_{T} > 1 GeV)", 100, -0.5, 249.5);
0189   chi2ndf = iBooker.book1D("chi2ndf", "PV tracks (p_{T} > 1 GeV) #chi^{2}/ndof", 100, 0., 20.);
0190   chi2prob = iBooker.book1D("chi2prob", "PV tracks (p_{T} > 1 GeV) #chi^{2} probability", 100, 0., 1.);
0191 
0192   dxy2 = iBooker.book1D("dxyzoom", "PV tracks (p_{T} > 1 GeV) d_{xy} (#mum)", DxyBin, DxyMin / 5., DxyMax / 5.);
0193 
0194   trackpt = pvMonitor::makeTH1IfLog(iBooker,
0195                                     true,
0196                                     false,
0197                                     "pt_track",
0198                                     "PV tracks p_{T};PV tracks p_{T} [GeV];#tracks",
0199                                     PtBin,
0200                                     log10(PtMin),
0201                                     log10(PtMax));
0202 
0203   phi_pt1 = iBooker.book1D("phi_pt1", "PV tracks (p_{T} > 1 GeV) #phi; PV tracks #phi;#tracks", PhiBin, PhiMin, PhiMax);
0204   eta_pt1 = iBooker.book1D("eta_pt1", "PV tracks (p_{T} > 1 GeV) #eta; PV tracks #eta;#tracks", EtaBin, EtaMin, EtaMax);
0205   phi_pt10 =
0206       iBooker.book1D("phi_pt10", "PV tracks (p_{T} > 10 GeV) #phi; PV tracks #phi;#tracks", PhiBin, PhiMin, PhiMax);
0207   eta_pt10 =
0208       iBooker.book1D("eta_pt10", "PV tracks (p_{T} > 10 GeV) #phi; PV tracks #eta;#tracks", EtaBin, EtaMin, EtaMax);
0209 
0210   // initialize and book the monitors;
0211   dxy_pt1.varname_ = "xy";
0212   dxy_pt1.pTcut_ = 1.f;
0213   dxy_pt1.bookIPMonitor(iBooker, conf_);
0214 
0215   dxy_pt10.varname_ = "xy";
0216   dxy_pt10.pTcut_ = 10.f;
0217   dxy_pt10.bookIPMonitor(iBooker, conf_);
0218 
0219   dz_pt1.varname_ = "z";
0220   dz_pt1.pTcut_ = 1.f;
0221   dz_pt1.bookIPMonitor(iBooker, conf_);
0222 
0223   dz_pt10.varname_ = "z";
0224   dz_pt10.pTcut_ = 10.f;
0225   dz_pt10.bookIPMonitor(iBooker, conf_);
0226 }
0227 
0228 void PrimaryVertexMonitor::IPMonitoring::bookIPMonitor(DQMStore::IBooker& iBooker, const edm::ParameterSet& config) {
0229   int VarBin = config.getParameter<int>(fmt::format("D{}Bin", varname_));
0230   double VarMin = config.getParameter<double>(fmt::format("D{}Min", varname_));
0231   double VarMax = config.getParameter<double>(fmt::format("D{}Max", varname_));
0232 
0233   PhiBin_ = config.getParameter<int>("PhiBin");
0234   PhiMin_ = config.getParameter<double>("PhiMin");
0235   PhiMax_ = config.getParameter<double>("PhiMax");
0236   int PhiBin2D = config.getParameter<int>("PhiBin2D");
0237 
0238   EtaBin_ = config.getParameter<int>("EtaBin");
0239   EtaMin_ = config.getParameter<double>("EtaMin");
0240   EtaMax_ = config.getParameter<double>("EtaMax");
0241   int EtaBin2D = config.getParameter<int>("EtaBin2D");
0242 
0243   PtBin_ = config.getParameter<int>("PtBin");
0244   PtMin_ = config.getParameter<double>("PtMin") * pTcut_;
0245   PtMax_ = config.getParameter<double>("PtMax") * pTcut_;
0246 
0247   // 1D variables
0248 
0249   IP_ = iBooker.book1D(fmt::format("d{}_pt{}", varname_, pTcut_),
0250                        fmt::format("PV tracks (p_{{T}} > {} GeV) d_{{{}}} (#mum)", pTcut_, varname_),
0251                        VarBin,
0252                        VarMin,
0253                        VarMax);
0254 
0255   IPErr_ = iBooker.book1D(fmt::format("d{}Err_pt{}", varname_, pTcut_),
0256                           fmt::format("PV tracks (p_{{T}} > {} GeV) d_{{{}}} error (#mum)", pTcut_, varname_),
0257                           100,
0258                           0.,
0259                           (varname_.find("xy") != std::string::npos) ? 2000. : 10000.);
0260 
0261   IPPull_ = iBooker.book1D(
0262       fmt::format("d{}Pull_pt{}", varname_, pTcut_),
0263       fmt::format("PV tracks (p_{{T}} > {} GeV) d_{{{}}}/#sigma_{{d_{{{}}}}}", pTcut_, varname_, varname_),
0264       100,
0265       -5.,
0266       5.);
0267 
0268   // IP profiles
0269 
0270   IPVsPhi_ = iBooker.bookProfile(fmt::format("d{}VsPhi_pt{}", varname_, pTcut_),
0271                                  fmt::format("PV tracks (p_{{T}} > {}) d_{{{}}} VS track #phi", pTcut_, varname_),
0272                                  PhiBin_,
0273                                  PhiMin_,
0274                                  PhiMax_,
0275                                  VarBin,
0276                                  VarMin,
0277                                  VarMax,
0278                                  "");
0279   IPVsPhi_->setAxisTitle("PV track (p_{T} > 1 GeV) #phi", 1);
0280   IPVsPhi_->setAxisTitle(fmt::format("PV tracks (p_{{T}} > {} GeV) d_{{{}}} (#mum)", pTcut_, varname_), 2);
0281 
0282   IPVsEta_ = iBooker.bookProfile(fmt::format("d{}VsEta_pt{}", varname_, pTcut_),
0283                                  fmt::format("PV tracks (p_{{T}} > {}) d_{{{}}} VS track #eta", pTcut_, varname_),
0284                                  EtaBin_,
0285                                  EtaMin_,
0286                                  EtaMax_,
0287                                  VarBin,
0288                                  VarMin,
0289                                  VarMax,
0290                                  "");
0291   IPVsEta_->setAxisTitle("PV track (p_{T} > 1 GeV) #eta", 1);
0292   IPVsEta_->setAxisTitle(fmt::format("PV tracks (p_{{T}} > {} GeV) d_{{{}}} (#mum)", pTcut_, varname_), 2);
0293 
0294   IPVsPt_ = pvMonitor::makeProfileIfLog(
0295       iBooker,
0296       true,  /* x-axis */
0297       false, /* y-axis */
0298       fmt::format("d{}VsPt_pt{}", varname_, pTcut_).c_str(),
0299       fmt::format("PV tracks (p_{{T}} > {}) d_{{{}}} VS track p_{{T}}", pTcut_, varname_).c_str(),
0300       PtBin_,
0301       log10(PtMin_),
0302       log10(PtMax_),
0303       VarMin,
0304       VarMax,
0305       "");
0306   IPVsPt_->setAxisTitle("PV track (p_{T} > 1 GeV) p_{T} [GeV]", 1);
0307   IPVsPt_->setAxisTitle(fmt::format("PV tracks (p_{{T}} > {} GeV) d_{{{}}} (#mum)", pTcut_, varname_), 2);
0308 
0309   // IP error profiles
0310 
0311   IPErrVsPhi_ =
0312       iBooker.bookProfile(fmt::format("d{}ErrVsPhi_pt{}", varname_, pTcut_),
0313                           fmt::format("PV tracks (p_{{T}} > {}) d_{{{}}} error VS track #phi", pTcut_, varname_),
0314                           PhiBin_,
0315                           PhiMin_,
0316                           PhiMax_,
0317                           VarBin,
0318                           0.,
0319                           (varname_.find("xy") != std::string::npos) ? 100. : 200.,
0320                           "");
0321   IPErrVsPhi_->setAxisTitle("PV track (p_{T} > 1 GeV) #phi", 1);
0322   IPErrVsPhi_->setAxisTitle(fmt::format("PV tracks (p_{{T}} > {} GeV) d_{{{}}} error (#mum)", pTcut_, varname_), 2);
0323 
0324   IPErrVsEta_ =
0325       iBooker.bookProfile(fmt::format("d{}ErrVsEta_pt{}", varname_, pTcut_),
0326                           fmt::format("PV tracks (p_{{T}} > {}) d_{{{}}} error VS track #eta", pTcut_, varname_),
0327                           EtaBin_,
0328                           EtaMin_,
0329                           EtaMax_,
0330                           VarBin,
0331                           0.,
0332                           (varname_.find("xy") != std::string::npos) ? 100. : 200.,
0333                           "");
0334   IPErrVsEta_->setAxisTitle("PV track (p_{T} > 1 GeV) #eta", 1);
0335   IPErrVsEta_->setAxisTitle(fmt::format("PV tracks (p_{{T}} > {} GeV) d_{{{}}} error (#mum)", pTcut_, varname_), 2);
0336 
0337   IPErrVsPt_ = pvMonitor::makeProfileIfLog(
0338       iBooker,
0339       true,  /* x-axis */
0340       false, /* y-axis */
0341       fmt::format("d{}ErrVsPt_pt{}", varname_, pTcut_).c_str(),
0342       fmt::format("PV tracks (p_{{T}} > {}) d_{{{}}} error VS track p_{{T}}", pTcut_, varname_).c_str(),
0343       PtBin_,
0344       log10(PtMin_),
0345       log10(PtMax_),
0346       VarMin,
0347       VarMax,
0348       "");
0349   IPErrVsPt_->setAxisTitle("PV track (p_{T} > 1 GeV) p_{T} [GeV]", 1);
0350   IPErrVsPt_->setAxisTitle(fmt::format("PV tracks (p_{{T}} > {} GeV) d_{{{}}} error (#mum)", pTcut_, varname_), 2);
0351 
0352   // 2D profiles
0353 
0354   IPVsEtaVsPhi_ = iBooker.bookProfile2D(
0355       fmt::format("d{}VsEtaVsPhi_pt{}", varname_, pTcut_),
0356       fmt::format("PV tracks (p_{{T}} > {}) d_{{{}}} VS track #eta VS track #phi", pTcut_, varname_),
0357       EtaBin2D,
0358       EtaMin_,
0359       EtaMax_,
0360       PhiBin2D,
0361       PhiMin_,
0362       PhiMax_,
0363       VarBin,
0364       VarMin,
0365       VarMax,
0366       "");
0367   IPVsEtaVsPhi_->setAxisTitle("PV track (p_{T} > 1 GeV) #eta", 1);
0368   IPVsEtaVsPhi_->setAxisTitle("PV track (p_{T} > 1 GeV) #phi", 2);
0369   IPVsEtaVsPhi_->setAxisTitle(fmt::format("PV tracks (p_{{T}} > {} GeV) d_{{{}}} (#mum)", pTcut_, varname_), 3);
0370 
0371   IPErrVsEtaVsPhi_ = iBooker.bookProfile2D(
0372       fmt::format("d{}ErrVsEtaVsPhi_pt{}", varname_, pTcut_),
0373       fmt::format("PV tracks (p_{{T}} > {}) d_{{{}}} error VS track #eta VS track #phi", pTcut_, varname_),
0374       EtaBin2D,
0375       EtaMin_,
0376       EtaMax_,
0377       PhiBin2D,
0378       PhiMin_,
0379       PhiMax_,
0380       VarBin,
0381       0.,
0382       (varname_.find("xy") != std::string::npos) ? 100. : 200.,
0383       "");
0384   IPErrVsEtaVsPhi_->setAxisTitle("PV track (p_{T} > 1 GeV) #eta", 1);
0385   IPErrVsEtaVsPhi_->setAxisTitle("PV track (p_{T} > 1 GeV) #phi", 2);
0386   IPErrVsEtaVsPhi_->setAxisTitle(fmt::format("PV tracks (p_{{T}} > {} GeV) d_{{{}}} error (#mum)", pTcut_, varname_),
0387                                  3);
0388 }
0389 
0390 void PrimaryVertexMonitor::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0391   Handle<reco::VertexCollection> recVtxs;
0392   iEvent.getByToken(vertexToken_, recVtxs);
0393 
0394   Handle<VertexScore> scores;
0395   iEvent.getByToken(scoreToken_, scores);
0396 
0397   edm::Handle<reco::BeamSpot> beamSpotHandle;
0398   iEvent.getByToken(beamspotToken_, beamSpotHandle);
0399 
0400   //
0401   // check for absent products and simply "return" in that case
0402   //
0403   if (recVtxs.isValid() == false || beamSpotHandle.isValid() == false) {
0404     edm::LogWarning("PrimaryVertexMonitor")
0405         << " Some products not available in the event: VertexCollection " << vertexInputTag_ << " " << recVtxs.isValid()
0406         << " BeamSpot " << beamSpotInputTag_ << " " << beamSpotHandle.isValid() << ". Skipping plots for this event";
0407     return;
0408   }
0409 
0410   // check upfront that refs to track are (likely) to be valid
0411   {
0412     bool ok = true;
0413     for (const auto& v : *recVtxs) {
0414       if (v.tracksSize() > 0) {
0415         const auto& ref = v.trackRefAt(0);
0416         if (ref.isNull() || !ref.isAvailable()) {
0417           if (!errorPrinted_)
0418             edm::LogWarning("PrimaryVertexMonitor")
0419                 << "Skipping vertex collection: " << vertexInputTag_
0420                 << " since likely the track collection the vertex has refs pointing to is missing (at least the first "
0421                    "TrackBaseRef is null or not available)";
0422           else
0423             errorPrinted_ = true;
0424           ok = false;
0425         }
0426       }
0427     }
0428     if (!ok)
0429       return;
0430   }
0431 
0432   BeamSpot beamSpot = *beamSpotHandle;
0433 
0434   nbvtx->Fill(recVtxs->size() * 1.);
0435   int ng = 0;
0436   for (auto const& vx : (*recVtxs))
0437     if (vx.isValid() && !vx.isFake() && vx.ndof() >= ndof_)
0438       ++ng;
0439   nbgvtx->Fill(ng * 1.);
0440 
0441   if (scores.isValid() && !(*scores).empty()) {
0442     auto pvScore = (*scores).get(0);
0443     score[1]->Fill(std::sqrt(pvScore));
0444     for (unsigned int i = 1; i < (*scores).size(); ++i)
0445       score[0]->Fill(std::sqrt((*scores).get(i)));
0446   }
0447 
0448   // fill PV tracks MEs (as now, for alignment)
0449   if (!recVtxs->empty()) {
0450     vertexPlots(recVtxs->front(), beamSpot, 1);
0451     pvTracksPlots(recVtxs->front());
0452 
0453     for (reco::VertexCollection::const_iterator v = recVtxs->begin() + 1; v != recVtxs->end(); ++v)
0454       vertexPlots(*v, beamSpot, 0);
0455   }
0456 
0457   // Beamline plots:
0458   bsX->Fill(beamSpot.x0());
0459   bsY->Fill(beamSpot.y0());
0460   bsZ->Fill(beamSpot.z0());
0461   bsSigmaZ->Fill(beamSpot.sigmaZ());
0462   bsDxdz->Fill(beamSpot.dxdz());
0463   bsDydz->Fill(beamSpot.dydz());
0464   bsBeamWidthX->Fill(beamSpot.BeamWidthX() * cmToUm);
0465   bsBeamWidthY->Fill(beamSpot.BeamWidthY() * cmToUm);
0466   bsType->Fill(beamSpot.type());
0467 }
0468 
0469 void PrimaryVertexMonitor::pvTracksPlots(const Vertex& v) {
0470   if (!v.isValid())
0471     return;
0472   if (v.isFake())
0473     return;
0474 
0475   if (v.tracksSize() == 0) {
0476     ntracks->Fill(0);
0477     return;
0478   }
0479 
0480   const math::XYZPoint myVertex(v.position().x(), v.position().y(), v.position().z());
0481 
0482   size_t nTracks = 0;
0483   float sumPT = 0.;
0484 
0485   for (reco::Vertex::trackRef_iterator t = v.tracks_begin(); t != v.tracks_end(); t++) {
0486     bool isHighPurity = (**t).quality(reco::TrackBase::highPurity);
0487     if (!isHighPurity && useHPfoAlignmentPlots_)
0488       continue;
0489 
0490     float pt = (**t).pt();
0491     trackpt->Fill(pt);
0492 
0493     if (pt < 1.)
0494       continue;
0495 
0496     nTracks++;
0497 
0498     float eta = (**t).eta();
0499     float phi = (**t).phi();
0500 
0501     float w = v.trackWeight(*t);
0502     float chi2NDF = (**t).normalizedChi2();
0503     float chi2Prob = TMath::Prob((**t).chi2(), (int)(**t).ndof());
0504     float Dxy = (**t).dxy(myVertex) * cmToUm;
0505     float Dz = (**t).dz(myVertex) * cmToUm;
0506     float DxyErr = (**t).dxyError() * cmToUm;
0507     float DzErr = (**t).dzError() * cmToUm;
0508 
0509     sumPT += pt * pt;
0510 
0511     // fill MEs
0512     phi_pt1->Fill(phi);
0513     eta_pt1->Fill(eta);
0514 
0515     weight->Fill(w);
0516     chi2ndf->Fill(chi2NDF);
0517     chi2prob->Fill(chi2Prob);
0518     dxy2->Fill(Dxy);
0519 
0520     // dxy pT>1
0521 
0522     dxy_pt1.IP_->Fill(Dxy);
0523     dxy_pt1.IPVsPhi_->Fill(phi, Dxy);
0524     dxy_pt1.IPVsEta_->Fill(eta, Dxy);
0525     dxy_pt1.IPVsPt_->Fill(pt, Dxy);
0526     dxy_pt1.IPVsEtaVsPhi_->Fill(eta, phi, Dxy);
0527 
0528     dxy_pt1.IPErr_->Fill(DxyErr);
0529     dxy_pt1.IPPull_->Fill(Dxy / DxyErr);
0530     dxy_pt1.IPErrVsPhi_->Fill(phi, DxyErr);
0531     dxy_pt1.IPErrVsEta_->Fill(eta, DxyErr);
0532     dxy_pt1.IPErrVsPt_->Fill(pt, DxyErr);
0533     dxy_pt1.IPErrVsEtaVsPhi_->Fill(eta, phi, DxyErr);
0534 
0535     // dz pT>1
0536 
0537     dz_pt1.IP_->Fill(Dz);
0538     dz_pt1.IPVsPhi_->Fill(phi, Dz);
0539     dz_pt1.IPVsEta_->Fill(eta, Dz);
0540     dz_pt1.IPVsPt_->Fill(pt, Dz);
0541     dz_pt1.IPVsEtaVsPhi_->Fill(eta, phi, Dz);
0542 
0543     dz_pt1.IPErr_->Fill(DzErr);
0544     dz_pt1.IPPull_->Fill(Dz / DzErr);
0545     dz_pt1.IPErrVsPhi_->Fill(phi, DzErr);
0546     dz_pt1.IPErrVsEta_->Fill(eta, DzErr);
0547     dz_pt1.IPErrVsPt_->Fill(pt, DxyErr);
0548     dz_pt1.IPErrVsEtaVsPhi_->Fill(eta, phi, DzErr);
0549 
0550     if (pt < 10.)
0551       continue;
0552 
0553     phi_pt10->Fill(phi);
0554     eta_pt10->Fill(eta);
0555 
0556     // dxy pT>10
0557 
0558     dxy_pt10.IP_->Fill(Dxy);
0559     dxy_pt10.IPVsPhi_->Fill(phi, Dxy);
0560     dxy_pt10.IPVsEta_->Fill(eta, Dxy);
0561     dxy_pt10.IPVsPt_->Fill(pt, Dxy);
0562     dxy_pt10.IPVsEtaVsPhi_->Fill(eta, phi, Dxy);
0563 
0564     dxy_pt10.IPErr_->Fill(DxyErr);
0565     dxy_pt10.IPPull_->Fill(Dxy / DxyErr);
0566     dxy_pt10.IPErrVsPhi_->Fill(phi, DxyErr);
0567     dxy_pt10.IPErrVsEta_->Fill(eta, DxyErr);
0568     dxy_pt10.IPErrVsPt_->Fill(pt, DxyErr);
0569     dxy_pt10.IPErrVsEtaVsPhi_->Fill(eta, phi, DxyErr);
0570 
0571     // dxz pT>10
0572 
0573     dz_pt10.IP_->Fill(Dz);
0574     dz_pt10.IPVsPhi_->Fill(phi, Dz);
0575     dz_pt10.IPVsEta_->Fill(eta, Dz);
0576     dz_pt10.IPVsPt_->Fill(pt, Dz);
0577     dz_pt10.IPVsEtaVsPhi_->Fill(eta, phi, Dz);
0578 
0579     dz_pt10.IPErr_->Fill(DzErr);
0580     dz_pt10.IPPull_->Fill(Dz / DzErr);
0581     dz_pt10.IPErrVsPhi_->Fill(phi, DzErr);
0582     dz_pt10.IPErrVsEta_->Fill(eta, DzErr);
0583     dz_pt10.IPErrVsPt_->Fill(pt, DxyErr);
0584     dz_pt10.IPErrVsEtaVsPhi_->Fill(eta, phi, DzErr);
0585   }
0586   ntracks->Fill(float(nTracks));
0587   sumpt->Fill(sumPT);
0588 }
0589 
0590 void PrimaryVertexMonitor::vertexPlots(const Vertex& v, const BeamSpot& beamSpot, int i) {
0591   if (i < 0 || i > 1)
0592     return;
0593   if (!v.isValid())
0594     type[i]->Fill(2.);
0595   else if (v.isFake())
0596     type[i]->Fill(1.);
0597   else
0598     type[i]->Fill(0.);
0599 
0600   if (v.isValid() && !v.isFake()) {
0601     float weight = 0;
0602     for (reco::Vertex::trackRef_iterator t = v.tracks_begin(); t != v.tracks_end(); t++)
0603       weight += v.trackWeight(*t);
0604     trksWeight[i]->Fill(weight);
0605     nbtksinvtx[i]->Fill(v.tracksSize());
0606     ntracksVsZ[i]->Fill(v.position().z() - beamSpot.z0(), v.tracksSize());
0607 
0608     vtxchi2[i]->Fill(v.chi2());
0609     vtxndf[i]->Fill(v.ndof());
0610     vtxprob[i]->Fill(ChiSquaredProbability(v.chi2(), v.ndof()));
0611 
0612     xrec[i]->Fill(v.position().x());
0613     yrec[i]->Fill(v.position().y());
0614     zrec[i]->Fill(v.position().z());
0615 
0616     float xb = beamSpot.x0() + beamSpot.dxdz() * (v.position().z() - beamSpot.z0());
0617     float yb = beamSpot.y0() + beamSpot.dydz() * (v.position().z() - beamSpot.z0());
0618     xDiff[i]->Fill((v.position().x() - xb) * cmToUm);
0619     yDiff[i]->Fill((v.position().y() - yb) * cmToUm);
0620 
0621     xerr[i]->Fill(v.xError() * cmToUm);
0622     yerr[i]->Fill(v.yError() * cmToUm);
0623     zerr[i]->Fill(v.zError() * cmToUm);
0624     xerrVsTrks[i]->Fill(weight, v.xError() * cmToUm);
0625     yerrVsTrks[i]->Fill(weight, v.yError() * cmToUm);
0626     zerrVsTrks[i]->Fill(weight, v.zError() * cmToUm);
0627 
0628     nans[i]->Fill(1., edm::isNotFinite(v.position().x()) * 1.);
0629     nans[i]->Fill(2., edm::isNotFinite(v.position().y()) * 1.);
0630     nans[i]->Fill(3., edm::isNotFinite(v.position().z()) * 1.);
0631 
0632     int index = 3;
0633     for (int k = 0; k != 3; k++) {
0634       for (int j = k; j != 3; j++) {
0635         index++;
0636         nans[i]->Fill(index * 1., edm::isNotFinite(v.covariance(k, j)) * 1.);
0637         // in addition, diagonal element must be positive
0638         if (j == k && v.covariance(k, j) < 0) {
0639           nans[i]->Fill(index * 1., 1.);
0640         }
0641       }
0642     }
0643   }
0644 }
0645 
0646 void PrimaryVertexMonitor::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0647   edm::ParameterSetDescription desc;
0648   desc.add<std::string>("TopFolderName", "OfflinePV");
0649   desc.add<std::string>("AlignmentLabel", "Alignment");
0650   desc.add<int>("ndof", 4);
0651   desc.add<bool>("useHPforAlignmentPlots", true);
0652   desc.add<InputTag>("vertexLabel", edm::InputTag("offlinePrimaryVertices"));
0653   desc.add<InputTag>("beamSpotLabel", edm::InputTag("offlineBeamSpot"));
0654   desc.add<double>("PUMax", 80.0);
0655   desc.add<double>("Xpos", 0.1);
0656   desc.add<double>("Ypos", 0.0);
0657   desc.add<int>("TkSizeBin", 100);
0658   desc.add<double>("TkSizeMin", -0.5);
0659   desc.add<double>("TkSizeMax", 499.5);
0660   desc.add<int>("DxyBin", 100);
0661   desc.add<double>("DxyMin", -5000.0);
0662   desc.add<double>("DxyMax", 5000.0);
0663   desc.add<int>("DzBin", 100);
0664   desc.add<double>("DzMin", -2000.0);
0665   desc.add<double>("DzMax", 2000.0);
0666   desc.add<int>("PhiBin", 32);
0667   desc.add<double>("PhiMin", -M_PI);
0668   desc.add<double>("PhiMax", M_PI);
0669   desc.add<int>("EtaBin", 26);
0670   desc.add<double>("EtaMin", 2.5);
0671   desc.add<double>("EtaMax", -2.5);
0672   desc.add<int>("PtBin", 49);
0673   desc.add<double>("PtMin", 1.);
0674   desc.add<double>("PtMax", 50.);
0675   desc.add<int>("PhiBin2D", 12);
0676   desc.add<int>("EtaBin2D", 8);
0677   descriptions.addWithDefaultLabel(desc);
0678 }
0679 
0680 //define this as a plug-in
0681 DEFINE_FWK_MODULE(PrimaryVertexMonitor);