File indexing completed on 2024-04-06 12:33:35
0001 #include "Validation/RecoVertex/interface/BSvsPVHistogramMaker.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0004 #include "FWCore/Framework/interface/Event.h"
0005 #include "FWCore/ServiceRegistry/interface/Service.h"
0006 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0007 #include "DataFormats/VertexReco/interface/Vertex.h"
0008 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0009 #include "TH1F.h"
0010 #include "TH2F.h"
0011 #include "TProfile.h"
0012
0013 BSvsPVHistogramMaker::BSvsPVHistogramMaker(edm::ConsumesCollector&& iC)
0014 : _currdir(nullptr),
0015 m_maxLS(100),
0016 useSlope_(true),
0017 _trueOnly(true),
0018 _runHisto(true),
0019 _runHistoProfile(true),
0020 _runHistoBXProfile(true),
0021 _runHistoBX2D(false),
0022 _histoParameters(),
0023 _rhm(iC) {}
0024
0025 BSvsPVHistogramMaker::BSvsPVHistogramMaker(const edm::ParameterSet& iConfig, edm::ConsumesCollector&& iC)
0026 : _currdir(nullptr),
0027 m_maxLS(iConfig.getParameter<unsigned int>("maxLSBeforeRebin")),
0028 useSlope_(iConfig.getParameter<bool>("useSlope")),
0029 _trueOnly(iConfig.getUntrackedParameter<bool>("trueOnly", true)),
0030 _runHisto(iConfig.getUntrackedParameter<bool>("runHisto", true)),
0031 _runHistoProfile(iConfig.getUntrackedParameter<bool>("runHistoProfile", true)),
0032 _runHistoBXProfile(iConfig.getUntrackedParameter<bool>("runHistoBXProfile", true)),
0033 _runHistoBX2D(iConfig.getUntrackedParameter<bool>("runHistoBX2D", false)),
0034 _histoParameters(iConfig.getUntrackedParameter<edm::ParameterSet>("histoParameters", edm::ParameterSet())),
0035 _rhm(iC) {}
0036
0037 BSvsPVHistogramMaker::~BSvsPVHistogramMaker() { delete _currdir; }
0038
0039 void BSvsPVHistogramMaker::book(const std::string dirname) {
0040 edm::Service<TFileService> tfserv;
0041 TFileDirectory* currdir = &(tfserv->tFileDirectory());
0042
0043 if (!dirname.empty()) {
0044 currdir = new TFileDirectory(tfserv->mkdir(dirname));
0045 _currdir = currdir;
0046 }
0047
0048 edm::LogInfo("HistogramBooking") << "Vertex histogram booking in directory " << dirname;
0049
0050 _hdeltax = currdir->make<TH1F>("deltax",
0051 "(PV-BS) X position",
0052 _histoParameters.getUntrackedParameter<unsigned int>("nBinX", 200),
0053 _histoParameters.getUntrackedParameter<double>("xMin", -1.),
0054 _histoParameters.getUntrackedParameter<double>("xMax", 1.));
0055 _hdeltax->GetXaxis()->SetTitle("#Delta(X) [cm]");
0056 _hdeltax->GetYaxis()->SetTitle("Vertices");
0057
0058 _hdeltay = currdir->make<TH1F>("deltay",
0059 "(PV-BS) Y position",
0060 _histoParameters.getUntrackedParameter<unsigned int>("nBinY", 200),
0061 _histoParameters.getUntrackedParameter<double>("yMin", -1.),
0062 _histoParameters.getUntrackedParameter<double>("yMax", 1.));
0063 _hdeltay->GetXaxis()->SetTitle("#Delta(Y) [cm]");
0064 _hdeltay->GetYaxis()->SetTitle("Vertices");
0065
0066 _hdeltaz = currdir->make<TH1F>("deltaz",
0067 "(PV-BS) Z position",
0068 _histoParameters.getUntrackedParameter<unsigned int>("nBinZ", 200),
0069 _histoParameters.getUntrackedParameter<double>("zMin", -20.),
0070 _histoParameters.getUntrackedParameter<double>("zMax", 20.));
0071 _hdeltaz->GetXaxis()->SetTitle("#Delta(Z) [cm]");
0072 _hdeltaz->GetYaxis()->SetTitle("Vertices");
0073
0074 _hdeltaxvsz = currdir->make<TProfile>("deltaxvsz",
0075 "(PV-BS) X position vs Z",
0076 _histoParameters.getUntrackedParameter<unsigned int>("nBinZProfile", 40),
0077 _histoParameters.getUntrackedParameter<double>("zMinProfile", -20.),
0078 _histoParameters.getUntrackedParameter<double>("zMaxProfile", 20.));
0079 _hdeltaxvsz->GetXaxis()->SetTitle("Z [cm]");
0080 _hdeltaxvsz->GetYaxis()->SetTitle("#Delta(X) [cm]");
0081
0082 _hdeltayvsz = currdir->make<TProfile>("deltayvsz",
0083 "(PV-BS) Y position vs Z",
0084 _histoParameters.getUntrackedParameter<unsigned int>("nBinZProfile", 40),
0085 _histoParameters.getUntrackedParameter<double>("zMinProfile", -20.),
0086 _histoParameters.getUntrackedParameter<double>("zMaxProfile", 20.));
0087 _hdeltayvsz->GetXaxis()->SetTitle("Z [cm]");
0088 _hdeltayvsz->GetYaxis()->SetTitle("#Delta(Y) [cm]");
0089
0090 if (_runHisto) {
0091 _hdeltaxrun = _rhm.makeTH1F("deltaxrun",
0092 "(PV-BS) X position",
0093 _histoParameters.getUntrackedParameter<unsigned int>("nBinX", 200),
0094 _histoParameters.getUntrackedParameter<double>("xMin", -1.),
0095 _histoParameters.getUntrackedParameter<double>("xMax", 1.));
0096
0097 _hdeltayrun = _rhm.makeTH1F("deltayrun",
0098 "(PV-BS) Y position",
0099 _histoParameters.getUntrackedParameter<unsigned int>("nBinY", 200),
0100 _histoParameters.getUntrackedParameter<double>("yMin", -1.),
0101 _histoParameters.getUntrackedParameter<double>("yMax", 1.));
0102
0103 _hdeltazrun = _rhm.makeTH1F("deltazrun",
0104 "(PV-BS) Z position",
0105 _histoParameters.getUntrackedParameter<unsigned int>("nBinZ", 200),
0106 _histoParameters.getUntrackedParameter<double>("zMin", -20.),
0107 _histoParameters.getUntrackedParameter<double>("zMax", 20.));
0108
0109 _hdeltaxvszrun = _rhm.makeTProfile("deltaxvszrun",
0110 "(PV-BS) X position vs Z",
0111 _histoParameters.getUntrackedParameter<unsigned int>("nBinZProfile", 40),
0112 _histoParameters.getUntrackedParameter<double>("zMinProfile", -20.),
0113 _histoParameters.getUntrackedParameter<double>("zMaxProfile", 20.));
0114
0115 _hdeltayvszrun = _rhm.makeTProfile("deltayvszrun",
0116 "(PV-BS) Y position vs Z",
0117 _histoParameters.getUntrackedParameter<unsigned int>("nBinZProfile", 40),
0118 _histoParameters.getUntrackedParameter<double>("zMinProfile", -20.),
0119 _histoParameters.getUntrackedParameter<double>("zMaxProfile", 20.));
0120
0121 if (_runHistoProfile) {
0122 _hdeltaxvsorbrun = _rhm.makeTProfile(
0123 "deltaxvsorbrun", "(PV-BS) X position vs orbit number", 4 * m_maxLS, 0.5, m_maxLS * 262144 + 0.5);
0124 _hdeltayvsorbrun = _rhm.makeTProfile(
0125 "deltayvsorbrun", "(PV-BS) Y position vs orbit number", 4 * m_maxLS, 0.5, m_maxLS * 262144 + 0.5);
0126 _hdeltazvsorbrun = _rhm.makeTProfile(
0127 "deltazvsorbrun", "(PV-BS) Z position vs orbit number", 4 * m_maxLS, 0.5, m_maxLS * 262144 + 0.5);
0128 }
0129 if (_runHistoBXProfile) {
0130 _hdeltaxvsbxrun = _rhm.makeTProfile("deltaxvsbxrun", "(PV-BS) X position vs BX number", 3564, -0.5, 3563.5);
0131 _hdeltayvsbxrun = _rhm.makeTProfile("deltayvsbxrun", "(PV-BS) Y position vs BX number", 3564, -0.5, 3563.5);
0132 _hdeltazvsbxrun = _rhm.makeTProfile("deltazvsbxrun", "(PV-BS) Z position vs BX number", 3564, -0.5, 3563.5);
0133 if (_runHistoBX2D) {
0134 _hdeltaxvsbx2drun = _rhm.makeTH2F("deltaxvsbx2drun",
0135 "(PV-BS) X position vs BX number",
0136 3564,
0137 -0.5,
0138 3563.5,
0139 _histoParameters.getUntrackedParameter<unsigned int>("nBinX", 200),
0140 _histoParameters.getUntrackedParameter<double>("xMin", -1.),
0141 _histoParameters.getUntrackedParameter<double>("xMax", 1.));
0142 _hdeltayvsbx2drun = _rhm.makeTH2F("deltayvsbx2drun",
0143 "(PV-BS) Y position vs BX number",
0144 3564,
0145 -0.5,
0146 3563.5,
0147 _histoParameters.getUntrackedParameter<unsigned int>("nBinY", 200),
0148 _histoParameters.getUntrackedParameter<double>("yMin", -1.),
0149 _histoParameters.getUntrackedParameter<double>("yMax", 1.));
0150 _hdeltazvsbx2drun = _rhm.makeTH2F("deltazvsbx2drun",
0151 "(PV-BS) Z position vs BX number",
0152 3564,
0153 -0.5,
0154 3563.5,
0155 _histoParameters.getUntrackedParameter<unsigned int>("nBinZ", 200),
0156 _histoParameters.getUntrackedParameter<double>("zMin", -20.),
0157 _histoParameters.getUntrackedParameter<double>("zMax", 20.));
0158 }
0159 }
0160 }
0161 }
0162
0163 void BSvsPVHistogramMaker::beginRun(const unsigned int nrun) {
0164 char runname[100];
0165 sprintf(runname, "run_%d", nrun);
0166
0167 TFileDirectory* currdir = _currdir;
0168 if (currdir == nullptr) {
0169 edm::Service<TFileService> tfserv;
0170 currdir = &(tfserv->tFileDirectory());
0171 }
0172
0173 _rhm.beginRun(nrun, *currdir);
0174
0175 if (_runHisto) {
0176 (*_hdeltaxrun)->GetXaxis()->SetTitle("#Delta(X) [cm]");
0177 (*_hdeltaxrun)->GetYaxis()->SetTitle("Vertices");
0178 (*_hdeltayrun)->GetXaxis()->SetTitle("#Delta(Y) [cm]");
0179 (*_hdeltayrun)->GetYaxis()->SetTitle("Vertices");
0180 (*_hdeltazrun)->GetXaxis()->SetTitle("#Delta(Z) [cm]");
0181 (*_hdeltazrun)->GetYaxis()->SetTitle("Vertices");
0182 (*_hdeltaxvszrun)->GetXaxis()->SetTitle("Z [cm]");
0183 (*_hdeltaxvszrun)->GetYaxis()->SetTitle("#Delta(X) [cm]");
0184 (*_hdeltayvszrun)->GetXaxis()->SetTitle("Z [cm]");
0185 (*_hdeltayvszrun)->GetYaxis()->SetTitle("#Delta(Y) [cm]");
0186
0187 if (_runHistoProfile) {
0188 (*_hdeltaxvsorbrun)->GetXaxis()->SetTitle("time [orbit#]");
0189 (*_hdeltaxvsorbrun)->GetYaxis()->SetTitle("#Delta(X) [cm]");
0190 (*_hdeltaxvsorbrun)->SetCanExtend(TH1::kAllAxes);
0191 (*_hdeltayvsorbrun)->GetXaxis()->SetTitle("time [orbit#]");
0192 (*_hdeltayvsorbrun)->GetYaxis()->SetTitle("#Delta(Y) [cm]");
0193 (*_hdeltayvsorbrun)->SetCanExtend(TH1::kAllAxes);
0194 (*_hdeltazvsorbrun)->GetXaxis()->SetTitle("time [orbit#]");
0195 (*_hdeltazvsorbrun)->GetYaxis()->SetTitle("#Delta(Z) [cm]");
0196 (*_hdeltazvsorbrun)->SetCanExtend(TH1::kAllAxes);
0197 }
0198 if (_runHistoBXProfile) {
0199 (*_hdeltaxvsbxrun)->GetXaxis()->SetTitle("BX");
0200 (*_hdeltaxvsbxrun)->GetYaxis()->SetTitle("#Delta(X) [cm]");
0201 (*_hdeltayvsbxrun)->GetXaxis()->SetTitle("BX");
0202 (*_hdeltayvsbxrun)->GetYaxis()->SetTitle("#Delta(Y) [cm]");
0203 (*_hdeltazvsbxrun)->GetXaxis()->SetTitle("BX");
0204 (*_hdeltazvsbxrun)->GetYaxis()->SetTitle("#Delta(Z) [cm]");
0205 if (_runHistoBX2D) {
0206 (*_hdeltaxvsbx2drun)->GetXaxis()->SetTitle("BX");
0207 (*_hdeltaxvsbx2drun)->GetYaxis()->SetTitle("#Delta(X) [cm]");
0208 (*_hdeltayvsbx2drun)->GetXaxis()->SetTitle("BX");
0209 (*_hdeltayvsbx2drun)->GetYaxis()->SetTitle("#Delta(Y) [cm]");
0210 (*_hdeltazvsbx2drun)->GetXaxis()->SetTitle("BX");
0211 (*_hdeltazvsbx2drun)->GetYaxis()->SetTitle("#Delta(Z) [cm]");
0212 }
0213 }
0214 }
0215 }
0216
0217 void BSvsPVHistogramMaker::fill(const unsigned int orbit,
0218 const int bx,
0219 const reco::VertexCollection& vertices,
0220 const reco::BeamSpot& bs) {
0221 for (reco::VertexCollection::const_iterator vtx = vertices.begin(); vtx != vertices.end(); ++vtx) {
0222 if (!(_trueOnly && vtx->isFake())) {
0223
0224
0225
0226
0227
0228 double deltax = vtx->x() - x(bs, vtx->z());
0229 double deltay = vtx->y() - y(bs, vtx->z());
0230 double deltaz = vtx->z() - bs.z0();
0231
0232 _hdeltax->Fill(deltax);
0233 _hdeltay->Fill(deltay);
0234 _hdeltaz->Fill(deltaz);
0235 _hdeltaxvsz->Fill(vtx->z(), deltax);
0236 _hdeltayvsz->Fill(vtx->z(), deltay);
0237
0238 if (_runHisto) {
0239 if (_hdeltaxrun && *_hdeltaxrun)
0240 (*_hdeltaxrun)->Fill(deltax);
0241 if (_hdeltayrun && *_hdeltayrun)
0242 (*_hdeltayrun)->Fill(deltay);
0243 if (_hdeltazrun && *_hdeltazrun)
0244 (*_hdeltazrun)->Fill(deltaz);
0245 if (_hdeltaxvszrun && *_hdeltaxvszrun)
0246 (*_hdeltaxvszrun)->Fill(vtx->z(), deltax);
0247 if (_hdeltayvszrun && *_hdeltayvszrun)
0248 (*_hdeltayvszrun)->Fill(vtx->z(), deltay);
0249 if (_runHistoProfile) {
0250 if (_hdeltaxvsorbrun && *_hdeltaxvsorbrun)
0251 (*_hdeltaxvsorbrun)->Fill(orbit, deltax);
0252 if (_hdeltayvsorbrun && *_hdeltayvsorbrun)
0253 (*_hdeltayvsorbrun)->Fill(orbit, deltay);
0254 if (_hdeltazvsorbrun && *_hdeltazvsorbrun)
0255 (*_hdeltazvsorbrun)->Fill(orbit, deltaz);
0256 }
0257 if (_runHistoBXProfile) {
0258 if (_hdeltaxvsbxrun && *_hdeltaxvsbxrun)
0259 (*_hdeltaxvsbxrun)->Fill(bx % 3564, deltax);
0260 if (_hdeltayvsbxrun && *_hdeltayvsbxrun)
0261 (*_hdeltayvsbxrun)->Fill(bx % 3564, deltay);
0262 if (_hdeltazvsbxrun && *_hdeltazvsbxrun)
0263 (*_hdeltazvsbxrun)->Fill(bx % 3564, deltaz);
0264 if (_runHistoBX2D) {
0265 if (_hdeltaxvsbx2drun && *_hdeltaxvsbx2drun)
0266 (*_hdeltaxvsbx2drun)->Fill(bx % 3564, deltax);
0267 if (_hdeltayvsbx2drun && *_hdeltayvsbx2drun)
0268 (*_hdeltayvsbx2drun)->Fill(bx % 3564, deltay);
0269 if (_hdeltazvsbx2drun && *_hdeltazvsbx2drun)
0270 (*_hdeltazvsbx2drun)->Fill(bx % 3564, deltaz);
0271 }
0272 }
0273 }
0274 }
0275 }
0276 }
0277
0278 void BSvsPVHistogramMaker::fill(const edm::Event& iEvent,
0279 const reco::VertexCollection& vertices,
0280 const reco::BeamSpot& bs) {
0281 fill(iEvent.orbitNumber(), iEvent.bunchCrossing(), vertices, bs);
0282 }
0283
0284 double BSvsPVHistogramMaker::x(const reco::BeamSpot& bs, const double z) const {
0285 double x = bs.x0();
0286
0287
0288 if (useSlope_)
0289 x += bs.dxdz() * (z - bs.z0());
0290
0291
0292
0293 return x;
0294 }
0295
0296 double BSvsPVHistogramMaker::y(const reco::BeamSpot& bs, const double z) const {
0297 double y = bs.y0();
0298
0299
0300 if (useSlope_)
0301 y += bs.dydz() * (z - bs.z0());
0302
0303
0304
0305 return y;
0306 }