Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:09:47

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