File indexing completed on 2024-04-06 12:07:05
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010 #include "DQM/DTMonitorClient/src/DTOccupancyTestML.h"
0011
0012 #include "FWCore/ServiceRegistry/interface/Service.h"
0013 #include "FWCore/Framework/interface/LuminosityBlock.h"
0014 #include "FWCore/Framework/interface/EventSetup.h"
0015 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0016 #include "Geometry/DTGeometry/interface/DTGeometry.h"
0017 #include "DQMServices/Core/interface/DQMStore.h"
0018 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0019 #include "FWCore/ParameterSet/interface/FileInPath.h"
0020
0021 #include "TMath.h"
0022
0023 #include "PhysicsTools/TensorFlow/interface/TensorFlow.h"
0024 #include <vector>
0025
0026 using namespace edm;
0027 using namespace std;
0028
0029 DTOccupancyTestML::DTOccupancyTestML(const edm::ParameterSet& ps)
0030 : muonGeomToken_(esConsumes<edm::Transition::BeginRun>()) {
0031 LogVerbatim("DTDQM|DTMonitorClient|DTOccupancyTestML") << "[DTOccupancyTestML]: Constructor";
0032
0033
0034
0035 lsCounter = 0;
0036
0037 writeRootFile = ps.getUntrackedParameter<bool>("writeRootFile", false);
0038 if (writeRootFile) {
0039 rootFile = new TFile("MLDTOccupancyTest.root", "RECREATE");
0040 ntuple = new TNtuple("OccupancyNtuple",
0041 "OccupancyNtuple",
0042 "ls:wh:st:se:lay1MeanCell:lay1RMS:lay2MeanCell:lay2RMS:lay3MeanCell:lay3RMS:lay4MeanCell:"
0043 "lay4RMS:lay5MeanCell:lay5RMS:lay6MeanCell:lay6RMS:lay7MeanCell:lay7RMS:lay8MeanCell:lay8RMS:"
0044 "lay9MeanCell:lay9RMS:lay10MeanCell:lay10RMS:lay11MeanCell:lay11RMS:lay12MeanCell:lay12RMS");
0045 }
0046
0047
0048 tpMode = ps.getUntrackedParameter<bool>("testPulseMode", false);
0049
0050 runOnAllHitsOccupancies = ps.getUntrackedParameter<bool>("runOnAllHitsOccupancies", true);
0051 runOnNoiseOccupancies = ps.getUntrackedParameter<bool>("runOnNoiseOccupancies", false);
0052 runOnInTimeOccupancies = ps.getUntrackedParameter<bool>("runOnInTimeOccupancies", false);
0053 nMinEvts = ps.getUntrackedParameter<int>("nEventsCert", 5000);
0054 nMinEvtsPC = ps.getUntrackedParameter<int>("nEventsMinPC", 2200);
0055 nZeroEvtsPC = ps.getUntrackedParameter<int>("nEventsZeroPC", 30);
0056
0057 bookingdone = false;
0058
0059
0060 nevents = 0;
0061 }
0062
0063 DTOccupancyTestML::~DTOccupancyTestML() {
0064 LogVerbatim("DTDQM|DTMonitorClient|MLDTOccupancyTest") << " destructor called" << endl;
0065 }
0066
0067 void DTOccupancyTestML::beginRun(const edm::Run& run, const EventSetup& context) {
0068 LogVerbatim("DTDQM|DTMonitorClient|DTOccupancyTestML") << "[DTOccupancyTestML]: BeginRun";
0069
0070
0071 muonGeom = &context.getData(muonGeomToken_);
0072 }
0073
0074 void DTOccupancyTestML::dqmEndLuminosityBlock(DQMStore::IBooker& ibooker,
0075 DQMStore::IGetter& igetter,
0076 edm::LuminosityBlock const& lumiSeg,
0077 edm::EventSetup const& context) {
0078 if (!bookingdone) {
0079
0080
0081 for (int wh = -2; wh <= 2; ++wh) {
0082 bookHistos(ibooker, wh, string("MLDTOccupancies"), "MLOccupancySummary");
0083 }
0084
0085 ibooker.setCurrentFolder(topFolder(true));
0086 string title = "Occupancy Summary";
0087 if (tpMode) {
0088 title = "Test Pulse Occupancy Summary";
0089 }
0090
0091 summaryHisto = ibooker.book2D("MLOccupancySummary", title.c_str(), 12, 1, 13, 5, -2, 3);
0092 summaryHisto->setAxisTitle("sector", 1);
0093 summaryHisto->setAxisTitle("wheel", 2);
0094
0095
0096 glbSummaryHisto = ibooker.book2D("MLOccupancyGlbSummary", title.c_str(), 12, 1, 13, 5, -2, 3);
0097 glbSummaryHisto->setAxisTitle("sector", 1);
0098 glbSummaryHisto->setAxisTitle("wheel", 2);
0099
0100
0101 if (runOnAllHitsOccupancies) {
0102 nameMonitoredHisto = "OccupancyAllHits_perCh";
0103 } else if (runOnNoiseOccupancies) {
0104 nameMonitoredHisto = "OccupancyNoise_perCh";
0105 } else if (runOnInTimeOccupancies) {
0106 nameMonitoredHisto = "OccupancyInTimeHits_perCh";
0107 } else {
0108 nameMonitoredHisto = "OccupancyAllHits_perCh";
0109 }
0110 }
0111 bookingdone = true;
0112
0113 LogVerbatim("DTDQM|DTMonitorClient|DTOccupancyTestML")
0114 << "[DTOccupancyTestML]: End of LS transition, performing the DQM client operation";
0115 lsCounter++;
0116
0117
0118 summaryHisto->Reset();
0119 glbSummaryHisto->Reset();
0120
0121 nChannelTotal = 0;
0122 nChannelDead = 0;
0123
0124
0125 vector<const DTChamber*> chambers = muonGeom->chambers();
0126
0127
0128 tensorflow::setLogging("3");
0129 edm::FileInPath modelFilePath("DQM/DTMonitorClient/data/occupancy_cnn_v1.pb");
0130 tensorflow::GraphDef* graphDef = tensorflow::loadGraphDef(modelFilePath.fullPath());
0131
0132
0133 tensorflow::Session* session = tensorflow::createSession(graphDef);
0134
0135 for (vector<const DTChamber*>::const_iterator chamber = chambers.begin(); chamber != chambers.end();
0136 ++chamber) {
0137 DTChamberId chId = (*chamber)->id();
0138
0139 MonitorElement* chamberOccupancyHisto = igetter.get(getMEName(nameMonitoredHisto, chId));
0140
0141
0142 if (chamberOccupancyHisto != nullptr) {
0143
0144 TH2F* histo = chamberOccupancyHisto->getTH2F();
0145
0146 float chamberPercentage = 1.;
0147 int result = runOccupancyTest(histo, chId, chamberPercentage, graphDef, session);
0148 int sector = chId.sector();
0149
0150 if (sector == 13) {
0151 sector = 4;
0152 float resultSect4 = wheelHistos[chId.wheel()]->getBinContent(sector, chId.station());
0153 if (resultSect4 > result) {
0154 result = (int)resultSect4;
0155 }
0156 } else if (sector == 14) {
0157 sector = 10;
0158 float resultSect10 = wheelHistos[chId.wheel()]->getBinContent(sector, chId.station());
0159 if (resultSect10 > result) {
0160 result = (int)resultSect10;
0161 }
0162 }
0163
0164
0165 if ((sector == 4 || sector == 10) && chId.station() == 4)
0166 chamberPercentage = chamberPercentage / 2.;
0167
0168 wheelHistos[chId.wheel()]->setBinContent(sector, chId.station(), result);
0169 if (result > summaryHisto->getBinContent(sector, chId.wheel() + 3)) {
0170 summaryHisto->setBinContent(sector, chId.wheel() + 3, result);
0171 }
0172 glbSummaryHisto->Fill(sector, chId.wheel(), chamberPercentage * 1. / 4.);
0173 } else {
0174 LogVerbatim("DTDQM|DTMonitorClient|DTOccupancyTestML")
0175 << "[DTOccupancyTestML] ME: " << getMEName(nameMonitoredHisto, chId) << " not found!" << endl;
0176 }
0177 }
0178
0179
0180 tensorflow::closeSession(session);
0181 delete graphDef;
0182
0183 string nEvtsName = "DT/EventInfo/Counters/nProcessedEventsDigi";
0184
0185 MonitorElement* meProcEvts = igetter.get(nEvtsName);
0186
0187 if (meProcEvts) {
0188 int nProcEvts = meProcEvts->getFloatValue();
0189 glbSummaryHisto->setEntries(nProcEvts < nMinEvts ? 10. : nProcEvts);
0190 summaryHisto->setEntries(nProcEvts < nMinEvts ? 10. : nProcEvts);
0191 } else {
0192 glbSummaryHisto->setEntries(nMinEvts + 1);
0193 summaryHisto->setEntries(nMinEvts + 1);
0194 LogVerbatim("DTDQM|DTMonitorClient|DTOccupancyTestML")
0195 << "[DTOccupancyTestML] ME: " << nEvtsName << " not found!" << endl;
0196 }
0197
0198
0199
0200
0201
0202 if (writeRootFile)
0203 ntuple->AutoSave("SaveSelf");
0204 }
0205
0206 void DTOccupancyTestML::dqmEndJob(DQMStore::IBooker& ibooker, DQMStore::IGetter& igetter) {
0207 LogVerbatim("DTDQM|DTMonitorClient|DTOccupancyTestML") << "[DTOccupancyTestML] endjob called!";
0208 if (writeRootFile) {
0209 rootFile->cd();
0210 ntuple->Write();
0211 rootFile->Close();
0212 }
0213 }
0214
0215
0216
0217 void DTOccupancyTestML::bookHistos(DQMStore::IBooker& ibooker, const int wheelId, string folder, string histoTag) {
0218
0219 stringstream wheel;
0220 wheel << wheelId;
0221
0222 ibooker.setCurrentFolder(topFolder(true));
0223
0224
0225 string histoName = histoTag + "_W" + wheel.str();
0226
0227 LogVerbatim("DTDQM|DTMonitorClient|DTOccupancyTestML")
0228 << "[DTOccupancyTestML]: booking wheel histo:" << histoName << " (tag " << histoTag
0229 << ") in: " << topFolder(true) + "Wheel" + wheel.str() + "/" + folder << endl;
0230
0231 string histoTitle = "Occupancy summary WHEEL: " + wheel.str();
0232 if (tpMode) {
0233 histoTitle = "TP Occupancy summary WHEEL: " + wheel.str();
0234 }
0235
0236 wheelHistos[wheelId] = ibooker.book2D(histoName, histoTitle, 12, 1, 13, 4, 1, 5);
0237 wheelHistos[wheelId]->setBinLabel(1, "MB1", 2);
0238 wheelHistos[wheelId]->setBinLabel(2, "MB2", 2);
0239 wheelHistos[wheelId]->setBinLabel(3, "MB3", 2);
0240 wheelHistos[wheelId]->setBinLabel(4, "MB4", 2);
0241 wheelHistos[wheelId]->setAxisTitle("sector", 1);
0242 }
0243
0244 string DTOccupancyTestML::getMEName(string histoTag, const DTChamberId& chId) {
0245 stringstream wheel;
0246 wheel << chId.wheel();
0247 stringstream station;
0248 station << chId.station();
0249 stringstream sector;
0250 sector << chId.sector();
0251
0252 string folderRoot =
0253 topFolder(false) + "Wheel" + wheel.str() + "/Sector" + sector.str() + "/Station" + station.str() + "/";
0254
0255
0256 string histoName = histoTag + "_W" + wheel.str() + "_St" + station.str() + "_Sec" + sector.str();
0257
0258 string histoname = folderRoot + histoName;
0259
0260 return histoname;
0261 }
0262
0263 int DTOccupancyTestML::getIntegral(TH2F* histo, int firstBinX, int lastBinX, int firstBinY, int lastBinY, bool doall) {
0264 int sum = 0;
0265 for (Int_t i = firstBinX; i < lastBinX + 1; i++) {
0266 for (Int_t j = firstBinY; j < lastBinY + 1; j++) {
0267 if (histo->GetBinContent(i, j) > 0) {
0268 if (!doall)
0269 return 1;
0270 sum += histo->GetBinContent(i, j);
0271 }
0272 }
0273 }
0274
0275 return sum;
0276 }
0277
0278 int DTOccupancyTestML::runOccupancyTest(TH2F* histo,
0279 const DTChamberId& chId,
0280 float& chamberPercentage,
0281 tensorflow::GraphDef* graphDef,
0282 tensorflow::Session* session) {
0283 LogTrace("DTDQM|DTMonitorClient|DTOccupancyTestML") << "--- Occupancy test for chamber: " << chId << endl;
0284
0285
0286 int totalLayers = 0;
0287 int badLayers = 0;
0288
0289
0290 for (int superLayer = 1; superLayer <= 3; superLayer++) {
0291 int binYlow = ((superLayer - 1) * 4) + 1;
0292
0293
0294 if (chId.station() == 4 && superLayer == 2)
0295 continue;
0296
0297
0298 for (int layer = 1; layer <= 4; layer++) {
0299 DTLayerId layID(chId, superLayer, layer);
0300 int firstWire = muonGeom->layer(layID)->specificTopology().firstChannel();
0301 int nWires = muonGeom->layer(layID)->specificTopology().channels();
0302 int binY = binYlow + (layer - 1);
0303 std::vector<float> layerOccupancy(nWires);
0304 int channelId = 0;
0305
0306
0307 for (int cell = firstWire; cell != (nWires + firstWire); cell++) {
0308 double cellOccupancy = histo->GetBinContent(cell, binY);
0309 layerOccupancy.at(channelId) = cellOccupancy;
0310 channelId++;
0311 }
0312
0313 int targetSize = 47;
0314 std::vector<float> reshapedLayerOccupancy = interpolateLayers(layerOccupancy, nWires, targetSize);
0315
0316
0317 float maxOccupancyInLayer = *std::max_element(reshapedLayerOccupancy.begin(), reshapedLayerOccupancy.end());
0318
0319 if (maxOccupancyInLayer != 0) {
0320 for (auto& val : reshapedLayerOccupancy)
0321 val /= maxOccupancyInLayer;
0322 }
0323
0324
0325 tensorflow::Tensor input(tensorflow::DT_FLOAT, {1, targetSize});
0326 for (int i = 0; i < targetSize; i++)
0327 input.matrix<float>()(0, i) = reshapedLayerOccupancy[i];
0328
0329 std::vector<tensorflow::Tensor> outputs;
0330 tensorflow::run(session, {{"input_cnn_input", input}}, {"output_cnn/Softmax"}, &outputs);
0331
0332 totalLayers++;
0333 bool isBad = outputs[0].matrix<float>()(0, 1) > 0.95;
0334 if (isBad)
0335 badLayers++;
0336 }
0337 }
0338
0339
0340 chamberPercentage = 1.0 - static_cast<float>(badLayers) / totalLayers;
0341
0342 if (badLayers > 8)
0343 return 3;
0344 if (badLayers > 4)
0345 return 2;
0346 if (badLayers > 0)
0347 return 1;
0348 return 0;
0349 }
0350
0351 std::vector<float> DTOccupancyTestML::interpolateLayers(std::vector<float> const& inputs, int size, int targetSize) {
0352
0353 int interpolationFloor = 0;
0354 float interpolationPoint = 0.;
0355 float step = static_cast<float>(size) / targetSize;
0356 std::vector<float> reshapedInput(targetSize);
0357
0358 int limitInLoop = inputs.size();
0359 limitInLoop = std::min(limitInLoop, targetSize) - 1;
0360 for (int i = 0; i < limitInLoop; i++) {
0361 interpolationFloor = static_cast<int>(std::floor(interpolationPoint));
0362
0363 reshapedInput.at(i) =
0364 (interpolationPoint - interpolationFloor) * (inputs[interpolationFloor + 1] - inputs[interpolationFloor]) +
0365 inputs[interpolationFloor];
0366 interpolationPoint = step + interpolationPoint;
0367 }
0368 return reshapedInput;
0369 }
0370
0371 string DTOccupancyTestML::topFolder(bool isBooking) const {
0372 if (tpMode)
0373 return string("DT/10-TestPulses/");
0374 if (isBooking)
0375 return string("DT/01-Digi/ML");
0376 return string("DT/01-Digi/");
0377 }