Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-05-20 22:25:04 

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 (!beamSpotHandle.isValid()) {
0404     edm::LogWarning("PrimaryVertexMonitor") << " Some products not available in the event: BeamSpot ("
0405                                             << beamSpotInputTag_ << "). Skipping plots for this event.";
0406     return;
0407   }
0408 
0409   if (!recVtxs.isValid()) {
0410     edm::LogWarning("PrimaryVertexMonitor") << " Some products not available in the event: VertexCollection ("
0411                                             << vertexInputTag_ << "). Skipping plots for this event.";
0412     return;
0413   }
0414 
0415   // check upfront that refs to track are (likely) to be valid
0416   {
0417     bool ok = true;
0418     for (const auto& v : *recVtxs) {
0419       if (v.tracksSize() > 0) {
0420         const auto& ref = v.trackRefAt(0);
0421         if (ref.isNull() || !ref.isAvailable()) {
0422           if (!errorPrinted_)
0423             edm::LogWarning("PrimaryVertexMonitor")
0424                 << "Skipping vertex collection: " << vertexInputTag_
0425                 << " since likely the track collection the vertex has refs pointing to is missing (at least the first "
0426                    "TrackBaseRef is null or not available)";
0427           else
0428             errorPrinted_ = true;
0429           ok = false;
0430         }
0431       }
0432     }
0433     if (!ok)
0434       return;
0435   }
0436 
0437   BeamSpot beamSpot = *beamSpotHandle;
0438 
0439   nbvtx->Fill(recVtxs->size() * 1.);
0440   int ng = 0;
0441   for (auto const& vx : (*recVtxs))
0442     if (vx.isValid() && !vx.isFake() && vx.ndof() >= ndof_)
0443       ++ng;
0444   nbgvtx->Fill(ng * 1.);
0445 
0446   if (scores.isValid() && !(*scores).empty()) {
0447     auto pvScore = (*scores).get(0);
0448     score[1]->Fill(std::sqrt(pvScore));
0449     for (unsigned int i = 1; i < (*scores).size(); ++i)
0450       score[0]->Fill(std::sqrt((*scores).get(i)));
0451   }
0452 
0453   // fill PV tracks MEs (as now, for alignment)
0454   if (!recVtxs->empty()) {
0455     vertexPlots(recVtxs->front(), beamSpot, 1);
0456     pvTracksPlots(recVtxs->front());
0457 
0458     for (reco::VertexCollection::const_iterator v = recVtxs->begin() + 1; v != recVtxs->end(); ++v)
0459       vertexPlots(*v, beamSpot, 0);
0460   }
0461 
0462   // Beamline plots:
0463   bsX->Fill(beamSpot.x0());
0464   bsY->Fill(beamSpot.y0());
0465   bsZ->Fill(beamSpot.z0());
0466   bsSigmaZ->Fill(beamSpot.sigmaZ());
0467   bsDxdz->Fill(beamSpot.dxdz());
0468   bsDydz->Fill(beamSpot.dydz());
0469   bsBeamWidthX->Fill(beamSpot.BeamWidthX() * cmToUm);
0470   bsBeamWidthY->Fill(beamSpot.BeamWidthY() * cmToUm);
0471   bsType->Fill(beamSpot.type());
0472 }
0473 
0474 void PrimaryVertexMonitor::pvTracksPlots(const Vertex& v) {
0475   if (!v.isValid())
0476     return;
0477   if (v.isFake())
0478     return;
0479 
0480   if (v.tracksSize() == 0) {
0481     ntracks->Fill(0);
0482     return;
0483   }
0484 
0485   const math::XYZPoint myVertex(v.position().x(), v.position().y(), v.position().z());
0486 
0487   size_t nTracks = 0;
0488   float sumPT = 0.;
0489 
0490   for (reco::Vertex::trackRef_iterator t = v.tracks_begin(); t != v.tracks_end(); t++) {
0491     bool isHighPurity = (**t).quality(reco::TrackBase::highPurity);
0492     if (!isHighPurity && useHPfoAlignmentPlots_)
0493       continue;
0494 
0495     float pt = (**t).pt();
0496     trackpt->Fill(pt);
0497 
0498     if (pt < 1.)
0499       continue;
0500 
0501     nTracks++;
0502 
0503     float eta = (**t).eta();
0504     float phi = (**t).phi();
0505 
0506     float w = v.trackWeight(*t);
0507     float chi2NDF = (**t).normalizedChi2();
0508     float chi2Prob = TMath::Prob((**t).chi2(), (int)(**t).ndof());
0509     float Dxy = (**t).dxy(myVertex) * cmToUm;
0510     float Dz = (**t).dz(myVertex) * cmToUm;
0511     float DxyErr = (**t).dxyError() * cmToUm;
0512     float DzErr = (**t).dzError() * cmToUm;
0513 
0514     sumPT += pt * pt;
0515 
0516     // fill MEs
0517     phi_pt1->Fill(phi);
0518     eta_pt1->Fill(eta);
0519 
0520     weight->Fill(w);
0521     chi2ndf->Fill(chi2NDF);
0522     chi2prob->Fill(chi2Prob);
0523     dxy2->Fill(Dxy);
0524 
0525     // dxy pT>1
0526 
0527     dxy_pt1.IP_->Fill(Dxy);
0528     dxy_pt1.IPVsPhi_->Fill(phi, Dxy);
0529     dxy_pt1.IPVsEta_->Fill(eta, Dxy);
0530     dxy_pt1.IPVsPt_->Fill(pt, 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.IPErrVsPt_->Fill(pt, DxyErr);
0538     dxy_pt1.IPErrVsEtaVsPhi_->Fill(eta, phi, DxyErr);
0539 
0540     // dz pT>1
0541 
0542     dz_pt1.IP_->Fill(Dz);
0543     dz_pt1.IPVsPhi_->Fill(phi, Dz);
0544     dz_pt1.IPVsEta_->Fill(eta, Dz);
0545     dz_pt1.IPVsPt_->Fill(pt, Dz);
0546     dz_pt1.IPVsEtaVsPhi_->Fill(eta, phi, Dz);
0547 
0548     dz_pt1.IPErr_->Fill(DzErr);
0549     dz_pt1.IPPull_->Fill(Dz / DzErr);
0550     dz_pt1.IPErrVsPhi_->Fill(phi, DzErr);
0551     dz_pt1.IPErrVsEta_->Fill(eta, DzErr);
0552     dz_pt1.IPErrVsPt_->Fill(pt, DzErr);
0553     dz_pt1.IPErrVsEtaVsPhi_->Fill(eta, phi, DzErr);
0554 
0555     if (pt < 10.)
0556       continue;
0557 
0558     phi_pt10->Fill(phi);
0559     eta_pt10->Fill(eta);
0560 
0561     // dxy pT>10
0562 
0563     dxy_pt10.IP_->Fill(Dxy);
0564     dxy_pt10.IPVsPhi_->Fill(phi, Dxy);
0565     dxy_pt10.IPVsEta_->Fill(eta, Dxy);
0566     dxy_pt10.IPVsPt_->Fill(pt, Dxy);
0567     dxy_pt10.IPVsEtaVsPhi_->Fill(eta, phi, Dxy);
0568 
0569     dxy_pt10.IPErr_->Fill(DxyErr);
0570     dxy_pt10.IPPull_->Fill(Dxy / DxyErr);
0571     dxy_pt10.IPErrVsPhi_->Fill(phi, DxyErr);
0572     dxy_pt10.IPErrVsEta_->Fill(eta, DxyErr);
0573     dxy_pt10.IPErrVsPt_->Fill(pt, DxyErr);
0574     dxy_pt10.IPErrVsEtaVsPhi_->Fill(eta, phi, DxyErr);
0575 
0576     // dxz pT>10
0577 
0578     dz_pt10.IP_->Fill(Dz);
0579     dz_pt10.IPVsPhi_->Fill(phi, Dz);
0580     dz_pt10.IPVsEta_->Fill(eta, Dz);
0581     dz_pt10.IPVsPt_->Fill(pt, Dz);
0582     dz_pt10.IPVsEtaVsPhi_->Fill(eta, phi, Dz);
0583 
0584     dz_pt10.IPErr_->Fill(DzErr);
0585     dz_pt10.IPPull_->Fill(Dz / DzErr);
0586     dz_pt10.IPErrVsPhi_->Fill(phi, DzErr);
0587     dz_pt10.IPErrVsEta_->Fill(eta, DzErr);
0588     dz_pt10.IPErrVsPt_->Fill(pt, DzErr);
0589     dz_pt10.IPErrVsEtaVsPhi_->Fill(eta, phi, DzErr);
0590   }
0591   ntracks->Fill(float(nTracks));
0592   sumpt->Fill(sumPT);
0593 }
0594 
0595 void PrimaryVertexMonitor::vertexPlots(const Vertex& v, const BeamSpot& beamSpot, int i) {
0596   if (i < 0 || i > 1)
0597     return;
0598   if (!v.isValid())
0599     type[i]->Fill(2.);
0600   else if (v.isFake())
0601     type[i]->Fill(1.);
0602   else
0603     type[i]->Fill(0.);
0604 
0605   if (v.isValid() && !v.isFake()) {
0606     float weight = 0;
0607     for (reco::Vertex::trackRef_iterator t = v.tracks_begin(); t != v.tracks_end(); t++)
0608       weight += v.trackWeight(*t);
0609     trksWeight[i]->Fill(weight);
0610     nbtksinvtx[i]->Fill(v.tracksSize());
0611     ntracksVsZ[i]->Fill(v.position().z() - beamSpot.z0(), v.tracksSize());
0612 
0613     vtxchi2[i]->Fill(v.chi2());
0614     vtxndf[i]->Fill(v.ndof());
0615     vtxprob[i]->Fill(ChiSquaredProbability(v.chi2(), v.ndof()));
0616 
0617     xrec[i]->Fill(v.position().x());
0618     yrec[i]->Fill(v.position().y());
0619     zrec[i]->Fill(v.position().z());
0620 
0621     float xb = beamSpot.x0() + beamSpot.dxdz() * (v.position().z() - beamSpot.z0());
0622     float yb = beamSpot.y0() + beamSpot.dydz() * (v.position().z() - beamSpot.z0());
0623     xDiff[i]->Fill((v.position().x() - xb) * cmToUm);
0624     yDiff[i]->Fill((v.position().y() - yb) * cmToUm);
0625 
0626     xerr[i]->Fill(v.xError() * cmToUm);
0627     yerr[i]->Fill(v.yError() * cmToUm);
0628     zerr[i]->Fill(v.zError() * cmToUm);
0629     xerrVsTrks[i]->Fill(weight, v.xError() * cmToUm);
0630     yerrVsTrks[i]->Fill(weight, v.yError() * cmToUm);
0631     zerrVsTrks[i]->Fill(weight, v.zError() * cmToUm);
0632 
0633     nans[i]->Fill(1., edm::isNotFinite(v.position().x()) * 1.);
0634     nans[i]->Fill(2., edm::isNotFinite(v.position().y()) * 1.);
0635     nans[i]->Fill(3., edm::isNotFinite(v.position().z()) * 1.);
0636 
0637     int index = 3;
0638     for (int k = 0; k != 3; k++) {
0639       for (int j = k; j != 3; j++) {
0640         index++;
0641         nans[i]->Fill(index * 1., edm::isNotFinite(v.covariance(k, j)) * 1.);
0642         // in addition, diagonal element must be positive
0643         if (j == k && v.covariance(k, j) < 0) {
0644           nans[i]->Fill(index * 1., 1.);
0645         }
0646       }
0647     }
0648   }
0649 }
0650 
0651 void PrimaryVertexMonitor::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0652   edm::ParameterSetDescription desc;
0653   desc.add<std::string>("TopFolderName", "OfflinePV");
0654   desc.add<std::string>("AlignmentLabel", "Alignment");
0655   desc.add<int>("ndof", 4);
0656   desc.add<bool>("useHPforAlignmentPlots", true);
0657   desc.add<InputTag>("vertexLabel", edm::InputTag("offlinePrimaryVertices"));
0658   desc.add<InputTag>("beamSpotLabel", edm::InputTag("offlineBeamSpot"));
0659   desc.add<double>("PUMax", 80.0);
0660   desc.add<double>("Xpos", 0.1);
0661   desc.add<double>("Ypos", 0.0);
0662   desc.add<int>("TkSizeBin", 100);
0663   desc.add<double>("TkSizeMin", -0.5);
0664   desc.add<double>("TkSizeMax", 499.5);
0665   desc.add<int>("DxyBin", 100);
0666   desc.add<double>("DxyMin", -5000.0);
0667   desc.add<double>("DxyMax", 5000.0);
0668   desc.add<int>("DzBin", 100);
0669   desc.add<double>("DzMin", -2000.0);
0670   desc.add<double>("DzMax", 2000.0);
0671   desc.add<int>("PhiBin", 32);
0672   desc.add<double>("PhiMin", -M_PI);
0673   desc.add<double>("PhiMax", M_PI);
0674   desc.add<int>("EtaBin", 26);
0675   desc.add<double>("EtaMin", 2.5);
0676   desc.add<double>("EtaMax", -2.5);
0677   desc.add<int>("PtBin", 49);
0678   desc.add<double>("PtMin", 1.);
0679   desc.add<double>("PtMax", 50.);
0680   desc.add<int>("PhiBin2D", 12);
0681   desc.add<int>("EtaBin2D", 8);
0682   descriptions.addWithDefaultLabel(desc);
0683 }
0684 
0685 //define this as a plug-in
0686 DEFINE_FWK_MODULE(PrimaryVertexMonitor);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!