DTMonitorClient/src/DTOccupancyTestML.cc

0001 /*
0002  *  See header file for a description of this class.
0003  *
0004  *  \author G. Cerminara - University and INFN Torino
0005  *
0006  *  threadsafe version (//-) oct/nov 2014 - WATWanAbdullah -ncpp-um-my
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   // Get the DQM service
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   // switch on the mode for running on test pulses (different top folder)
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   // Event counter
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   // Get the geometry
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     // Book the summary histos
0080     //   - one summary per wheel
0081     for (int wh = -2; wh <= 2; ++wh) {  // loop over wheels
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     //   - global summary with alarms
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     //   - global summary with percentages
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     // assign the name of the input histogram
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 {  // default is AllHits histo
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   // Reset the global summary
0118   summaryHisto->Reset();
0119   glbSummaryHisto->Reset();
0120
0121   nChannelTotal = 0;
0122   nChannelDead = 0;
0123
0124   // Get all the DT chambers
0125   vector<const DTChamber*> chambers = muonGeom->chambers();
0126
0127   // Load graph
0128   edm::FileInPath modelFilePath("DQM/DTMonitorClient/data/occupancy_cnn_v1.pb");
0129   tensorflow::GraphDef* graphDef = tensorflow::loadGraphDef(modelFilePath.fullPath());
0130
0131   // Create session
0132   tensorflow::Session* session = tensorflow::createSession(graphDef);
0133
0134   for (vector<const DTChamber*>::const_iterator chamber = chambers.begin(); chamber != chambers.end();
0135        ++chamber) {  // Loop over all chambers
0136     DTChamberId chId = (*chamber)->id();
0137
0138     MonitorElement* chamberOccupancyHisto = igetter.get(getMEName(nameMonitoredHisto, chId));
0139
0140     // Run the tests on the plot for the various granularities
0141     if (chamberOccupancyHisto != nullptr) {
0142       // Get the 2D histo
0143       TH2F* histo = chamberOccupancyHisto->getTH2F();
0144
0145       float chamberPercentage = 1.;
0146       int result = runOccupancyTest(histo, chId, chamberPercentage, graphDef, session);
0147       int sector = chId.sector();
0148
0149       if (sector == 13) {
0150         sector = 4;
0151         float resultSect4 = wheelHistos[chId.wheel()]->getBinContent(sector, chId.station());
0152         if (resultSect4 > result) {
0153           result = (int)resultSect4;
0154         }
0155       } else if (sector == 14) {
0156         sector = 10;
0157         float resultSect10 = wheelHistos[chId.wheel()]->getBinContent(sector, chId.station());
0158         if (resultSect10 > result) {
0159           result = (int)resultSect10;
0160         }
0161       }
0162
0163       // the 2 MB4 of Sect 4 and 10 count as half a chamber
0164       if ((sector == 4 || sector == 10) && chId.station() == 4)
0165         chamberPercentage = chamberPercentage / 2.;
0166
0167       wheelHistos[chId.wheel()]->setBinContent(sector, chId.station(), result);
0168       if (result > summaryHisto->getBinContent(sector, chId.wheel() + 3)) {
0169         summaryHisto->setBinContent(sector, chId.wheel() + 3, result);
0170       }
0171       glbSummaryHisto->Fill(sector, chId.wheel(), chamberPercentage * 1. / 4.);
0172     } else {
0173       LogVerbatim("DTDQM|DTMonitorClient|DTOccupancyTestML")
0174           << "[DTOccupancyTestML] ME: " << getMEName(nameMonitoredHisto, chId) << " not found!" << endl;
0175     }
0176   }
0177
0178   // Clean up neural network graph
0179   tensorflow::closeSession(session);
0180   delete graphDef;
0181
0182   string nEvtsName = "DT/EventInfo/Counters/nProcessedEventsDigi";
0183
0184   MonitorElement* meProcEvts = igetter.get(nEvtsName);
0185
0186   if (meProcEvts) {
0187     int nProcEvts = meProcEvts->getFloatValue();
0188     glbSummaryHisto->setEntries(nProcEvts < nMinEvts ? 10. : nProcEvts);
0189     summaryHisto->setEntries(nProcEvts < nMinEvts ? 10. : nProcEvts);
0190   } else {
0191     glbSummaryHisto->setEntries(nMinEvts + 1);
0192     summaryHisto->setEntries(nMinEvts + 1);
0193     LogVerbatim("DTDQM|DTMonitorClient|DTOccupancyTestML")
0194         << "[DTOccupancyTestML] ME: " << nEvtsName << " not found!" << endl;
0195   }
0196
0197   // Fill the global summary
0198   // Check for entire sectors off and report them on the global summary
0199   //FIXME: TODO
0200
0201   if (writeRootFile)
0202     ntuple->AutoSave("SaveSelf");
0203 }
0204
0205 void DTOccupancyTestML::dqmEndJob(DQMStore::IBooker& ibooker, DQMStore::IGetter& igetter) {
0206   LogVerbatim("DTDQM|DTMonitorClient|DTOccupancyTestML") << "[DTOccupancyTestML] endjob called!";
0207   if (writeRootFile) {
0208     rootFile->cd();
0209     ntuple->Write();
0210     rootFile->Close();
0211   }
0212 }
0213
0214 // --------------------------------------------------
0215
0216 void DTOccupancyTestML::bookHistos(DQMStore::IBooker& ibooker, const int wheelId, string folder, string histoTag) {
0217   // Set the current folder
0218   stringstream wheel;
0219   wheel << wheelId;
0220
0221   ibooker.setCurrentFolder(topFolder(true));
0222
0223   // build the histo name
0224   string histoName = histoTag + "_W" + wheel.str();
0225
0226   LogVerbatim("DTDQM|DTMonitorClient|DTOccupancyTestML")
0227       << "[DTOccupancyTestML]: booking wheel histo:" << histoName << " (tag " << histoTag
0228       << ") in: " << topFolder(true) + "Wheel" + wheel.str() + "/" + folder << endl;
0229
0230   string histoTitle = "Occupancy summary WHEEL: " + wheel.str();
0231   if (tpMode) {
0232     histoTitle = "TP Occupancy summary WHEEL: " + wheel.str();
0233   }
0234
0235   wheelHistos[wheelId] = ibooker.book2D(histoName, histoTitle, 12, 1, 13, 4, 1, 5);
0236   wheelHistos[wheelId]->setBinLabel(1, "MB1", 2);
0237   wheelHistos[wheelId]->setBinLabel(2, "MB2", 2);
0238   wheelHistos[wheelId]->setBinLabel(3, "MB3", 2);
0239   wheelHistos[wheelId]->setBinLabel(4, "MB4", 2);
0240   wheelHistos[wheelId]->setAxisTitle("sector", 1);
0241 }
0242
0243 string DTOccupancyTestML::getMEName(string histoTag, const DTChamberId& chId) {
0244   stringstream wheel;
0245   wheel << chId.wheel();
0246   stringstream station;
0247   station << chId.station();
0248   stringstream sector;
0249   sector << chId.sector();
0250
0251   string folderRoot =
0252       topFolder(false) + "Wheel" + wheel.str() + "/Sector" + sector.str() + "/Station" + station.str() + "/";
0253
0254   // build the histo name
0255   string histoName = histoTag + "_W" + wheel.str() + "_St" + station.str() + "_Sec" + sector.str();
0256
0257   string histoname = folderRoot + histoName;
0258
0259   return histoname;
0260 }
0261
0262 int DTOccupancyTestML::getIntegral(TH2F* histo, int firstBinX, int lastBinX, int firstBinY, int lastBinY, bool doall) {
0263   int sum = 0;
0264   for (Int_t i = firstBinX; i < lastBinX + 1; i++) {
0265     for (Int_t j = firstBinY; j < lastBinY + 1; j++) {
0266       if (histo->GetBinContent(i, j) > 0) {
0267         if (!doall)
0268           return 1;
0269         sum += histo->GetBinContent(i, j);
0270       }
0271     }
0272   }
0273
0274   return sum;
0275 }
0276
0277 int DTOccupancyTestML::runOccupancyTest(TH2F* histo,
0278                                         const DTChamberId& chId,
0279                                         float& chamberPercentage,
0280                                         tensorflow::GraphDef* graphDef,
0281                                         tensorflow::Session* session) {
0282   LogTrace("DTDQM|DTMonitorClient|DTOccupancyTestML") << "--- Occupancy test for chamber: " << chId << endl;
0283
0284   // Initialize counters
0285   int totalLayers = 0;
0286   int badLayers = 0;
0287
0288   // Loop over the super layers
0289   for (int superLayer = 1; superLayer <= 3; superLayer++) {
0290     int binYlow = ((superLayer - 1) * 4) + 1;
0291
0292     // Skip for non-existent super layers
0293     if (chId.station() == 4 && superLayer == 2)
0294       continue;
0295
0296     // Loop over layers
0297     for (int layer = 1; layer <= 4; layer++) {
0298       DTLayerId layID(chId, superLayer, layer);
0299       int firstWire = muonGeom->layer(layID)->specificTopology().firstChannel();
0300       int nWires = muonGeom->layer(layID)->specificTopology().channels();
0301       int binY = binYlow + (layer - 1);
0302       std::vector<float> layerOccupancy(nWires);
0303       int channelId = 0;
0304
0305       // Loop over cells within  a layer
0306       for (int cell = firstWire; cell != (nWires + firstWire); cell++) {
0307         double cellOccupancy = histo->GetBinContent(cell, binY);
0308         layerOccupancy.at(channelId) = cellOccupancy;
0309         channelId++;
0310       }
0311
0312       int targetSize = 47;
0313       std::vector<float> reshapedLayerOccupancy = interpolateLayers(layerOccupancy, nWires, targetSize);
0314
0315       // Scale occupancy
0316       float maxOccupancyInLayer = *std::max_element(reshapedLayerOccupancy.begin(), reshapedLayerOccupancy.end());
0317
0318       if (maxOccupancyInLayer != 0) {
0319         for (auto& val : reshapedLayerOccupancy)
0320           val /= maxOccupancyInLayer;
0321       }
0322
0323       // Define input
0324       tensorflow::Tensor input(tensorflow::DT_FLOAT, {1, targetSize});
0325       for (int i = 0; i < targetSize; i++)
0326         input.matrix<float>()(0, i) = reshapedLayerOccupancy[i];
0327
0328       std::vector<tensorflow::Tensor> outputs;
0329       tensorflow::run(session, {{"input_cnn_input", input}}, {"output_cnn/Softmax"}, &outputs);
0330
0331       totalLayers++;
0332       bool isBad = outputs[0].matrix<float>()(0, 1) > 0.95;
0333       if (isBad)
0334         badLayers++;
0335     }
0336   }
0337
0338   // Calculate a fraction of good layers
0339   chamberPercentage = 1.0 - static_cast<float>(badLayers) / totalLayers;
0340
0341   if (badLayers > 8)
0342     return 3;  // 3 SL
0343   if (badLayers > 4)
0344     return 2;  // 2 SL
0345   if (badLayers > 0)
0346     return 1;  // 1 SL
0347   return 0;
0348 }
0349
0350 std::vector<float> DTOccupancyTestML::interpolateLayers(std::vector<float> const& inputs, int size, int targetSize) {
0351   // Reshape layers with linear interpolation
0352   int interpolationFloor = 0;
0353   float interpolationPoint = 0.;
0354   float step = static_cast<float>(size) / targetSize;
0355   std::vector<float> reshapedInput(targetSize);
0356
0357   int limitInLoop = inputs.size();
0358   limitInLoop = std::min(limitInLoop, targetSize) - 1;
0359   for (int i = 0; i < limitInLoop; i++) {
0360     interpolationFloor = static_cast<int>(std::floor(interpolationPoint));
0361     // Interpolating here
0362     reshapedInput.at(i) =
0363         (interpolationPoint - interpolationFloor) * (inputs[interpolationFloor + 1] - inputs[interpolationFloor]) +
0364         inputs[interpolationFloor];
0365     interpolationPoint = step + interpolationPoint;
0366   }
0367   return reshapedInput;
0368 }
0369
0370 string DTOccupancyTestML::topFolder(bool isBooking) const {
0371   if (tpMode)
0372     return string("DT/10-TestPulses/");
0373   if (isBooking)
0374     return string("DT/01-Digi/ML");
0375   return string("DT/01-Digi/");
0376 }