File indexing completed on 2022-05-14 00:30:14
0001
0002
0003
0004
0005
0006
0007 #include "FWCore/Framework/interface/MakerMacros.h"
0008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0009 #include "FWCore/Framework/interface/ESHandle.h"
0010 #include "FWCore/ServiceRegistry/interface/Service.h"
0011 #include "DQM/BeamMonitor/plugins/OnlineBeamMonitor.h"
0012 #include "DataFormats/Common/interface/Handle.h"
0013 #include "DataFormats/Common/interface/View.h"
0014 #include "FWCore/Framework/interface/ConsumesCollector.h"
0015 #include "FWCore/Framework/interface/EventSetup.h"
0016 #include "FWCore/Framework/interface/LuminosityBlock.h"
0017 #include "FWCore/Framework/interface/Run.h"
0018 #include "RecoVertex/BeamSpotProducer/interface/BeamFitter.h"
0019 #include "RecoVertex/BeamSpotProducer/interface/PVFitter.h"
0020 #include <memory>
0021
0022 #include <numeric>
0023
0024 using namespace std;
0025 using namespace edm;
0026 using namespace reco;
0027
0028
0029 OnlineBeamMonitor::OnlineBeamMonitor(const ParameterSet& ps)
0030 : monitorName_(ps.getUntrackedParameter<string>("MonitorName")),
0031 bsTransientToken_(esConsumes<edm::Transition::BeginLuminosityBlock>()),
0032 bsHLTToken_(esConsumes<edm::Transition::BeginLuminosityBlock>()),
0033 bsLegacyToken_(esConsumes<edm::Transition::BeginLuminosityBlock>()),
0034 numberOfValuesToSave_(0),
0035 appendRunTxt_(ps.getUntrackedParameter<bool>("AppendRunToFileName")),
0036 writeDIPTxt_(ps.getUntrackedParameter<bool>("WriteDIPAscii")),
0037 outputDIPTxt_(ps.getUntrackedParameter<std::string>("DIPFileName")) {
0038 if (!monitorName_.empty())
0039 monitorName_ = monitorName_ + "/";
0040
0041 processedLumis_.clear();
0042
0043 varNamesV_.push_back("x");
0044 varNamesV_.push_back("y");
0045 varNamesV_.push_back("z");
0046 varNamesV_.push_back("sigmaX");
0047 varNamesV_.push_back("sigmaY");
0048 varNamesV_.push_back("sigmaZ");
0049
0050
0051
0052
0053
0054
0055
0056
0057
0058 histoByCategoryNames_.insert(pair<string, string>("lumi", "Lumibased BeamSpotHLT"));
0059 histoByCategoryNames_.insert(pair<string, string>("lumi", "Lumibased BeamSpotLegacy"));
0060 histoByCategoryNames_.insert(pair<string, string>("lumi", "Lumibased BeamSpotTransient"));
0061
0062 for (const auto& itV : varNamesV_) {
0063 for (const auto& itM : histoByCategoryNames_) {
0064 histosMap_[itV][itM.first][itM.second] = nullptr;
0065 }
0066 }
0067 }
0068
0069 void OnlineBeamMonitor::fillDescriptions(edm::ConfigurationDescriptions& iDesc) {
0070 edm::ParameterSetDescription ps;
0071 ps.addUntracked<std::string>("MonitorName", "YourSubsystemName");
0072 ps.addUntracked<bool>("AppendRunToFileName", false);
0073 ps.addUntracked<bool>("WriteDIPAscii", true);
0074 ps.addUntracked<std::string>("DIPFileName", "BeamFitResultsForDIP.txt");
0075
0076 iDesc.addWithDefaultLabel(ps);
0077 }
0078
0079
0080 void OnlineBeamMonitor::bookHistograms(DQMStore::IBooker& ibooker,
0081 edm::Run const& iRun,
0082 edm::EventSetup const& iSetup) {
0083 string name;
0084 string title;
0085 int firstLumi = 1;
0086 int lastLumi = 3000;
0087 for (auto& itM : histosMap_) {
0088
0089
0090 for (auto& itMM : itM.second) {
0091 if (itMM.first != "run") {
0092 for (auto& itMMM : itMM.second) {
0093 name = string("h") + itM.first + itMMM.first;
0094 title = itM.first + "_{0} " + itMMM.first;
0095 if (itMM.first == "lumi") {
0096 ibooker.setCurrentFolder(monitorName_ + "Debug");
0097 itMMM.second = ibooker.book1D(name, title, lastLumi - firstLumi + 1, firstLumi - 0.5, lastLumi + 0.5);
0098 itMMM.second->setEfficiencyFlag();
0099 } else {
0100 LogInfo("OnlineBeamMonitorClient") << "Unrecognized category " << itMM.first;
0101 }
0102 if (itMMM.second != nullptr) {
0103 if (itMMM.first.find('-') != string::npos) {
0104 itMMM.second->setAxisTitle(string("#Delta ") + itM.first + "_{0} (cm)", 2);
0105 } else {
0106 itMMM.second->setAxisTitle(itM.first + "_{0} (cm)", 2);
0107 }
0108 itMMM.second->setAxisTitle("Lumisection", 1);
0109 }
0110 }
0111 }
0112 }
0113 }
0114
0115
0116 ibooker.setCurrentFolder(monitorName_ + "Validation");
0117
0118 bsChoice_ = ibooker.bookProfile("bsChoice",
0119 "BS Choice (+1): HLT - (-1): Legacy - (-10): Fake BS - (0): No Transient ",
0120 lastLumi - firstLumi + 1,
0121 firstLumi - 0.5,
0122 lastLumi + 0.5,
0123 100,
0124 -10,
0125 1,
0126 "");
0127 bsChoice_->setAxisTitle("Lumisection", 1);
0128 bsChoice_->setAxisTitle("Choice", 2);
0129 }
0130
0131
0132 std::shared_ptr<onlinebeammonitor::NoCache> OnlineBeamMonitor::globalBeginLuminosityBlock(
0133 const LuminosityBlock& iLumi, const EventSetup& iSetup) const {
0134
0135
0136 processedLumis_.push_back(iLumi.id().luminosityBlock());
0137
0138 ESHandle<BeamSpotOnlineObjects> bsHLTHandle;
0139 ESHandle<BeamSpotOnlineObjects> bsLegacyHandle;
0140 ESHandle<BeamSpotObjects> bsTransientHandle;
0141
0142
0143 std::string startTimeStamp_ = "0";
0144 std::string startTimeStampHLT_ = "0";
0145 std::string startTimeStampLegacy_ = "0";
0146 std::string stopTimeStamp_ = "0";
0147 std::string stopTimeStampHLT_ = "0";
0148 std::string stopTimeStampLegacy_ = "0";
0149 std::string lumiRange_ = "0 - 0";
0150 std::string lumiRangeHLT_ = "0 - 0";
0151 std::string lumiRangeLegacy_ = "0 - 0";
0152
0153 if (auto bsHLTHandle = iSetup.getHandle(bsHLTToken_)) {
0154 auto const& spotDB = *bsHLTHandle;
0155
0156
0157 startTimeStampHLT_ = spotDB.startTime();
0158 stopTimeStampHLT_ = spotDB.endTime();
0159 lumiRangeHLT_ = spotDB.lumiRange();
0160
0161
0162 BeamSpot::Point apoint(spotDB.x(), spotDB.y(), spotDB.z());
0163
0164 BeamSpot::CovarianceMatrix matrix;
0165 for (int i = 0; i < 7; ++i) {
0166 for (int j = 0; j < 7; ++j) {
0167 matrix(i, j) = spotDB.covariance(i, j);
0168 }
0169 }
0170
0171 beamSpotsMap_["HLT"] = BeamSpot(apoint, spotDB.sigmaZ(), spotDB.dxdz(), spotDB.dydz(), spotDB.beamWidthX(), matrix);
0172
0173 BeamSpot* aSpot = &(beamSpotsMap_["HLT"]);
0174
0175 aSpot->setBeamWidthY(spotDB.beamWidthY());
0176 aSpot->setEmittanceX(spotDB.emittanceX());
0177 aSpot->setEmittanceY(spotDB.emittanceY());
0178 aSpot->setbetaStar(spotDB.betaStar());
0179
0180 if (spotDB.beamType() == 2) {
0181 aSpot->setType(reco::BeamSpot::Tracker);
0182 } else {
0183 aSpot->setType(reco::BeamSpot::Fake);
0184 }
0185
0186
0187 } else {
0188 LogInfo("OnlineBeamMonitor") << "Database BeamSpot is not valid at lumi: " << iLumi.id().luminosityBlock();
0189 }
0190 if (auto bsLegacyHandle = iSetup.getHandle(bsLegacyToken_)) {
0191 auto const& spotDB = *bsLegacyHandle;
0192
0193
0194 BeamSpot::Point apoint(spotDB.x(), spotDB.y(), spotDB.z());
0195
0196
0197 startTimeStampLegacy_ = spotDB.startTime();
0198 stopTimeStampLegacy_ = spotDB.endTime();
0199 lumiRangeLegacy_ = spotDB.lumiRange();
0200
0201 BeamSpot::CovarianceMatrix matrix;
0202 for (int i = 0; i < 7; ++i) {
0203 for (int j = 0; j < 7; ++j) {
0204 matrix(i, j) = spotDB.covariance(i, j);
0205 }
0206 }
0207
0208 beamSpotsMap_["Legacy"] =
0209 BeamSpot(apoint, spotDB.sigmaZ(), spotDB.dxdz(), spotDB.dydz(), spotDB.beamWidthX(), matrix);
0210
0211 BeamSpot* aSpot = &(beamSpotsMap_["Legacy"]);
0212
0213 aSpot->setBeamWidthY(spotDB.beamWidthY());
0214 aSpot->setEmittanceX(spotDB.emittanceX());
0215 aSpot->setEmittanceY(spotDB.emittanceY());
0216 aSpot->setbetaStar(spotDB.betaStar());
0217
0218 if (spotDB.beamType() == 2) {
0219 aSpot->setType(reco::BeamSpot::Tracker);
0220 } else {
0221 aSpot->setType(reco::BeamSpot::Fake);
0222 }
0223
0224
0225 } else {
0226 LogInfo("OnlineBeamMonitor") << "Database BeamSpot is not valid at lumi: " << iLumi.id().luminosityBlock();
0227 }
0228 if (auto bsTransientHandle = iSetup.getHandle(bsTransientToken_)) {
0229 auto const& spotDB = *bsTransientHandle;
0230
0231
0232
0233 BeamSpot::Point apoint(spotDB.x(), spotDB.y(), spotDB.z());
0234
0235 BeamSpot::CovarianceMatrix matrix;
0236 for (int i = 0; i < 7; ++i) {
0237 for (int j = 0; j < 7; ++j) {
0238 matrix(i, j) = spotDB.covariance(i, j);
0239 }
0240 }
0241
0242 beamSpotsMap_["Transient"] =
0243 BeamSpot(apoint, spotDB.sigmaZ(), spotDB.dxdz(), spotDB.dydz(), spotDB.beamWidthX(), matrix);
0244
0245 BeamSpot* aSpot = &(beamSpotsMap_["Transient"]);
0246
0247 aSpot->setBeamWidthY(spotDB.beamWidthY());
0248 aSpot->setEmittanceX(spotDB.emittanceX());
0249 aSpot->setEmittanceY(spotDB.emittanceY());
0250 aSpot->setbetaStar(spotDB.betaStar());
0251 if (spotDB.beamType() == 2) {
0252 aSpot->setType(reco::BeamSpot::Tracker);
0253 } else {
0254 aSpot->setType(reco::BeamSpot::Fake);
0255 }
0256
0257 if (writeDIPTxt_) {
0258 std::ofstream outFile;
0259
0260 std::string tmpname = outputDIPTxt_;
0261 int frun = iLumi.getRun().run();
0262
0263 char index[15];
0264 if (appendRunTxt_ && writeDIPTxt_) {
0265 sprintf(index, "%s%i", "_Run", frun);
0266 tmpname.insert(outputDIPTxt_.length() - 4, index);
0267 }
0268
0269
0270 if (beamSpotsMap_.find("Transient") != beamSpotsMap_.end()) {
0271 if (beamSpotsMap_.find("HLT") != beamSpotsMap_.end() &&
0272 beamSpotsMap_["Transient"].x0() == beamSpotsMap_["HLT"].x0()) {
0273
0274 startTimeStamp_ = startTimeStampHLT_;
0275 stopTimeStamp_ = stopTimeStampHLT_;
0276 lumiRange_ = lumiRangeHLT_;
0277
0278 } else if (beamSpotsMap_.find("Legacy") != beamSpotsMap_.end() &&
0279 beamSpotsMap_["Transient"].x0() == beamSpotsMap_["Legacy"].x0()) {
0280
0281 startTimeStamp_ = startTimeStampLegacy_;
0282 stopTimeStamp_ = stopTimeStampLegacy_;
0283 lumiRange_ = lumiRangeLegacy_;
0284 }
0285 }
0286
0287 outFile.open(tmpname.c_str());
0288
0289 outFile << "Runnumber " << frun << " bx " << 0 << std::endl;
0290 outFile << "BeginTimeOfFit " << startTimeStamp_ << " " << 0 << std::endl;
0291 outFile << "EndTimeOfFit " << stopTimeStamp_ << " " << 0 << std::endl;
0292
0293 outFile << "LumiRange " << lumiRange_ << std::endl;
0294 outFile << "Type " << aSpot->type() << std::endl;
0295 outFile << "X0 " << aSpot->x0() << std::endl;
0296 outFile << "Y0 " << aSpot->y0() << std::endl;
0297 outFile << "Z0 " << aSpot->z0() << std::endl;
0298 outFile << "sigmaZ0 " << aSpot->sigmaZ() << std::endl;
0299 outFile << "dxdz " << aSpot->dxdz() << std::endl;
0300 outFile << "dydz " << aSpot->dydz() << std::endl;
0301 outFile << "BeamWidthX " << aSpot->BeamWidthX() << std::endl;
0302 outFile << "BeamWidthY " << aSpot->BeamWidthY() << std::endl;
0303 for (int i = 0; i < 6; ++i) {
0304 outFile << "Cov(" << i << ",j) ";
0305 for (int j = 0; j < 7; ++j) {
0306 outFile << aSpot->covariance(i, j) << " ";
0307 }
0308 outFile << std::endl;
0309 }
0310 outFile << "Cov(6,j) 0 0 0 0 0 0 " << aSpot->covariance(6, 6) << std::endl;
0311 outFile << "EmittanceX " << aSpot->emittanceX() << std::endl;
0312 outFile << "EmittanceY " << aSpot->emittanceY() << std::endl;
0313 outFile << "BetaStar " << aSpot->betaStar() << std::endl;
0314
0315 outFile.close();
0316 }
0317
0318
0319 } else {
0320 LogInfo("OnlineBeamMonitor") << "Database BeamSpot is not valid at lumi: " << iLumi.id().luminosityBlock();
0321 }
0322 return nullptr;
0323 }
0324
0325
0326 void OnlineBeamMonitor::globalEndLuminosityBlock(const LuminosityBlock& iLumi, const EventSetup& iSetup) {
0327
0328 if (beamSpotsMap_.find("Transient") != beamSpotsMap_.end()) {
0329 if (beamSpotsMap_.find("HLT") != beamSpotsMap_.end() &&
0330 beamSpotsMap_["Transient"].x0() == beamSpotsMap_["HLT"].x0()) {
0331 bsChoice_->Fill(iLumi.id().luminosityBlock(), 1);
0332 bsChoice_->setBinError(iLumi.id().luminosityBlock(), 0.05);
0333 } else if (beamSpotsMap_.find("Legacy") != beamSpotsMap_.end() &&
0334 beamSpotsMap_["Transient"].x0() == beamSpotsMap_["Legacy"].x0()) {
0335 bsChoice_->Fill(iLumi.id().luminosityBlock(), -1);
0336 bsChoice_->setBinError(iLumi.id().luminosityBlock(), 0.05);
0337 } else {
0338 bsChoice_->Fill(iLumi.id().luminosityBlock(), -10);
0339 bsChoice_->setBinError(iLumi.id().luminosityBlock(), 0.05);
0340 }
0341 } else {
0342 bsChoice_->Fill(iLumi.id().luminosityBlock(), 0);
0343 bsChoice_->setBinError(iLumi.id().luminosityBlock(), 0.05);
0344 }
0345
0346
0347 map<std::string, pair<double, double> > resultsMap;
0348 vector<pair<double, double> > vertexResults;
0349 MonitorElement* histo = nullptr;
0350 for (const auto& itV : varNamesV_) {
0351 resultsMap.clear();
0352 for (const auto& itBS : beamSpotsMap_) {
0353 if (itBS.second.type() == BeamSpot::Tracker) {
0354 if (itV == "x") {
0355 resultsMap[itBS.first] = pair<double, double>(itBS.second.x0(), itBS.second.x0Error());
0356 } else if (itV == "y") {
0357 resultsMap[itBS.first] = pair<double, double>(itBS.second.y0(), itBS.second.y0Error());
0358 } else if (itV == "z") {
0359 resultsMap[itBS.first] = pair<double, double>(itBS.second.z0(), itBS.second.z0Error());
0360 } else if (itV == "sigmaX") {
0361 resultsMap[itBS.first] = pair<double, double>(itBS.second.BeamWidthX(), itBS.second.BeamWidthXError());
0362 } else if (itV == "sigmaY") {
0363 resultsMap[itBS.first] = pair<double, double>(itBS.second.BeamWidthY(), itBS.second.BeamWidthYError());
0364 } else if (itV == "sigmaZ") {
0365 resultsMap[itBS.first] = pair<double, double>(itBS.second.sigmaZ(), itBS.second.sigmaZ0Error());
0366 } else {
0367 LogInfo("OnlineBeamMonitor") << "The histosMap_ has been built with the name " << itV
0368 << " that I can't recognize!";
0369 }
0370 }
0371 }
0372
0373 for (const auto& itM : histoByCategoryNames_) {
0374 if ((histo = histosMap_[itV][itM.first][itM.second]) == nullptr)
0375 continue;
0376 if (itM.second == "Lumibased BeamSpotHLT") {
0377 if (resultsMap.find("HLT") != resultsMap.end()) {
0378 histo->setBinContent(iLumi.id().luminosityBlock(), resultsMap["HLT"].first);
0379 histo->setBinError(iLumi.id().luminosityBlock(), resultsMap["HLT"].second);
0380 }
0381 } else if (itM.second == "Lumibased BeamSpotLegacy") {
0382 if (resultsMap.find("Legacy") != resultsMap.end()) {
0383 histo->setBinContent(iLumi.id().luminosityBlock(), resultsMap["Legacy"].first);
0384 histo->setBinError(iLumi.id().luminosityBlock(), resultsMap["Legacy"].second);
0385 }
0386 } else if (itM.second == "Lumibased BeamSpotTransient") {
0387 if (resultsMap.find("Transient") != resultsMap.end()) {
0388 histo->setBinContent(iLumi.id().luminosityBlock(), resultsMap["Transient"].first);
0389 histo->setBinError(iLumi.id().luminosityBlock(), resultsMap["Transient"].second);
0390 }
0391 } else {
0392 LogInfo("OnlineBeamMonitor") << "The histosMap_ have a histogram named " << itM.second
0393 << " that I can't recognize in this loop!";
0394 }
0395 }
0396 }
0397 }
0398
0399 void OnlineBeamMonitor::dqmEndRun(edm::Run const&, edm::EventSetup const&) {
0400 if (processedLumis_.empty()) {
0401 return;
0402 }
0403
0404 const double bigNumber = 1000000.;
0405 std::sort(processedLumis_.begin(), processedLumis_.end());
0406 int firstLumi = *processedLumis_.begin();
0407 int lastLumi = *(--processedLumis_.end());
0408 bsChoice_->getTH1()->GetXaxis()->SetRangeUser(firstLumi - 0.5, lastLumi + 0.5);
0409 for (auto& itH : histosMap_) {
0410 for (auto& itHH : itH.second) {
0411 double min = bigNumber;
0412 double max = -bigNumber;
0413 if (itHH.first != "run") {
0414 for (auto& itHHH : itHH.second) {
0415 if (itHHH.second != nullptr) {
0416 for (int bin = 1; bin <= itHHH.second->getTH1()->GetNbinsX(); bin++) {
0417 if (itHHH.second->getTH1()->GetBinError(bin) != 0 || itHHH.second->getTH1()->GetBinContent(bin) != 0) {
0418 if (itHHH.first == "Lumibased BeamSpotHLT" || itHHH.first == "Lumibased BeamSpotLegacy" ||
0419 itHHH.first == "Lumibased BeamSpotTransient") {
0420 if (min > itHHH.second->getTH1()->GetBinContent(bin)) {
0421 min = itHHH.second->getTH1()->GetBinContent(bin);
0422 }
0423 if (max < itHHH.second->getTH1()->GetBinContent(bin)) {
0424 max = itHHH.second->getTH1()->GetBinContent(bin);
0425 }
0426 } else {
0427 LogInfo("OnlineBeamMonitorClient") << "The histosMap_ have a histogram named " << itHHH.first
0428 << " that I can't recognize in this loop!";
0429 }
0430 }
0431 }
0432 }
0433 }
0434 for (auto& itHHH : itHH.second) {
0435 if (itHHH.second != nullptr) {
0436 if (itHHH.first == "Lumibased BeamSpotHLT" || itHHH.first == "Lumibased BeamSpotLegacy" ||
0437 itHHH.first == "Lumibased BeamSpotTransient") {
0438 if ((max == -bigNumber && min == bigNumber) || max - min == 0) {
0439 itHHH.second->getTH1()->SetMinimum(itHHH.second->getTH1()->GetMinimum() - 0.01);
0440 itHHH.second->getTH1()->SetMaximum(itHHH.second->getTH1()->GetMaximum() + 0.01);
0441 } else {
0442 itHHH.second->getTH1()->SetMinimum(min - 0.1 * (max - min));
0443 itHHH.second->getTH1()->SetMaximum(max + 0.1 * (max - min));
0444 }
0445 } else {
0446 LogInfo("OnlineBeamMonitorClient")
0447 << "The histosMap_ have a histogram named " << itHHH.first << " that I can't recognize in this loop!";
0448 }
0449 itHHH.second->getTH1()->GetXaxis()->SetRangeUser(firstLumi - 0.5, lastLumi + 0.5);
0450 }
0451 }
0452 }
0453 }
0454 }
0455 }
0456
0457 DEFINE_FWK_MODULE(OnlineBeamMonitor);