CMS 3D CMS Logo

DTOccupancyTestML.cc
Go to the documentation of this file.
1 /*
2  * See header file for a description of this class.
3  *
4  * \author G. Cerminara - University and INFN Torino
5  *
6  * threadsafe version (//-) oct/nov 2014 - WATWanAbdullah -ncpp-um-my
7  *
8  */
9 
11 
20 
21 #include "TMath.h"
22 
24 #include <vector>
25 
26 using namespace edm;
27 using namespace std;
28 
30  : muonGeomToken_(esConsumes<edm::Transition::BeginRun>()) {
31  LogVerbatim("DTDQM|DTMonitorClient|DTOccupancyTestML") << "[DTOccupancyTestML]: Constructor";
32 
33  // Get the DQM service
34 
35  lsCounter = 0;
36 
37  writeRootFile = ps.getUntrackedParameter<bool>("writeRootFile", false);
38  if (writeRootFile) {
39  rootFile = new TFile("MLDTOccupancyTest.root", "RECREATE");
40  ntuple = new TNtuple("OccupancyNtuple",
41  "OccupancyNtuple",
42  "ls:wh:st:se:lay1MeanCell:lay1RMS:lay2MeanCell:lay2RMS:lay3MeanCell:lay3RMS:lay4MeanCell:"
43  "lay4RMS:lay5MeanCell:lay5RMS:lay6MeanCell:lay6RMS:lay7MeanCell:lay7RMS:lay8MeanCell:lay8RMS:"
44  "lay9MeanCell:lay9RMS:lay10MeanCell:lay10RMS:lay11MeanCell:lay11RMS:lay12MeanCell:lay12RMS");
45  }
46 
47  // switch on the mode for running on test pulses (different top folder)
48  tpMode = ps.getUntrackedParameter<bool>("testPulseMode", false);
49 
50  runOnAllHitsOccupancies = ps.getUntrackedParameter<bool>("runOnAllHitsOccupancies", true);
51  runOnNoiseOccupancies = ps.getUntrackedParameter<bool>("runOnNoiseOccupancies", false);
52  runOnInTimeOccupancies = ps.getUntrackedParameter<bool>("runOnInTimeOccupancies", false);
53  nMinEvts = ps.getUntrackedParameter<int>("nEventsCert", 5000);
54  nMinEvtsPC = ps.getUntrackedParameter<int>("nEventsMinPC", 2200);
55  nZeroEvtsPC = ps.getUntrackedParameter<int>("nEventsZeroPC", 30);
56 
57  bookingdone = false;
58 
59  // Event counter
60  nevents = 0;
61 }
62 
64  LogVerbatim("DTDQM|DTMonitorClient|MLDTOccupancyTest") << " destructor called" << endl;
65 }
66 
68  LogVerbatim("DTDQM|DTMonitorClient|DTOccupancyTestML") << "[DTOccupancyTestML]: BeginRun";
69 
70  // Get the geometry
71  muonGeom = &context.getData(muonGeomToken_);
72 }
73 
75  DQMStore::IGetter& igetter,
76  edm::LuminosityBlock const& lumiSeg,
77  edm::EventSetup const& context) {
78  if (!bookingdone) {
79  // Book the summary histos
80  // - one summary per wheel
81  for (int wh = -2; wh <= 2; ++wh) { // loop over wheels
82  bookHistos(ibooker, wh, string("MLDTOccupancies"), "MLOccupancySummary");
83  }
84 
85  ibooker.setCurrentFolder(topFolder(true));
86  string title = "Occupancy Summary";
87  if (tpMode) {
88  title = "Test Pulse Occupancy Summary";
89  }
90  // - global summary with alarms
91  summaryHisto = ibooker.book2D("MLOccupancySummary", title.c_str(), 12, 1, 13, 5, -2, 3);
92  summaryHisto->setAxisTitle("sector", 1);
93  summaryHisto->setAxisTitle("wheel", 2);
94 
95  // - global summary with percentages
96  glbSummaryHisto = ibooker.book2D("MLOccupancyGlbSummary", title.c_str(), 12, 1, 13, 5, -2, 3);
97  glbSummaryHisto->setAxisTitle("sector", 1);
98  glbSummaryHisto->setAxisTitle("wheel", 2);
99 
100  // assign the name of the input histogram
102  nameMonitoredHisto = "OccupancyAllHits_perCh";
103  } else if (runOnNoiseOccupancies) {
104  nameMonitoredHisto = "OccupancyNoise_perCh";
105  } else if (runOnInTimeOccupancies) {
106  nameMonitoredHisto = "OccupancyInTimeHits_perCh";
107  } else { // default is AllHits histo
108  nameMonitoredHisto = "OccupancyAllHits_perCh";
109  }
110  }
111  bookingdone = true;
112 
113  LogVerbatim("DTDQM|DTMonitorClient|DTOccupancyTestML")
114  << "[DTOccupancyTestML]: End of LS transition, performing the DQM client operation";
115  lsCounter++;
116 
117  // Reset the global summary
118  summaryHisto->Reset();
120 
121  nChannelTotal = 0;
122  nChannelDead = 0;
123 
124  // Get all the DT chambers
125  vector<const DTChamber*> chambers = muonGeom->chambers();
126 
127  // Load graph
129  edm::FileInPath modelFilePath("DQM/DTMonitorClient/data/occupancy_cnn_v1.pb");
130  tensorflow::GraphDef* graphDef = tensorflow::loadGraphDef(modelFilePath.fullPath());
131 
132  // Create session
133  tensorflow::Session* session = tensorflow::createSession(graphDef);
134 
135  for (vector<const DTChamber*>::const_iterator chamber = chambers.begin(); chamber != chambers.end();
136  ++chamber) { // Loop over all chambers
137  DTChamberId chId = (*chamber)->id();
138 
139  MonitorElement* chamberOccupancyHisto = igetter.get(getMEName(nameMonitoredHisto, chId));
140 
141  // Run the tests on the plot for the various granularities
142  if (chamberOccupancyHisto != nullptr) {
143  // Get the 2D histo
144  TH2F* histo = chamberOccupancyHisto->getTH2F();
145 
146  float chamberPercentage = 1.;
147  int result = runOccupancyTest(histo, chId, chamberPercentage, graphDef, session);
148  int sector = chId.sector();
149 
150  if (sector == 13) {
151  sector = 4;
152  float resultSect4 = wheelHistos[chId.wheel()]->getBinContent(sector, chId.station());
153  if (resultSect4 > result) {
154  result = (int)resultSect4;
155  }
156  } else if (sector == 14) {
157  sector = 10;
158  float resultSect10 = wheelHistos[chId.wheel()]->getBinContent(sector, chId.station());
159  if (resultSect10 > result) {
160  result = (int)resultSect10;
161  }
162  }
163 
164  // the 2 MB4 of Sect 4 and 10 count as half a chamber
165  if ((sector == 4 || sector == 10) && chId.station() == 4)
166  chamberPercentage = chamberPercentage / 2.;
167 
168  wheelHistos[chId.wheel()]->setBinContent(sector, chId.station(), result);
169  if (result > summaryHisto->getBinContent(sector, chId.wheel() + 3)) {
171  }
172  glbSummaryHisto->Fill(sector, chId.wheel(), chamberPercentage * 1. / 4.);
173  } else {
174  LogVerbatim("DTDQM|DTMonitorClient|DTOccupancyTestML")
175  << "[DTOccupancyTestML] ME: " << getMEName(nameMonitoredHisto, chId) << " not found!" << endl;
176  }
177  }
178 
179  // Clean up neural network graph
180  tensorflow::closeSession(session);
181  delete graphDef;
182 
183  string nEvtsName = "DT/EventInfo/Counters/nProcessedEventsDigi";
184 
185  MonitorElement* meProcEvts = igetter.get(nEvtsName);
186 
187  if (meProcEvts) {
188  int nProcEvts = meProcEvts->getFloatValue();
189  glbSummaryHisto->setEntries(nProcEvts < nMinEvts ? 10. : nProcEvts);
190  summaryHisto->setEntries(nProcEvts < nMinEvts ? 10. : nProcEvts);
191  } else {
194  LogVerbatim("DTDQM|DTMonitorClient|DTOccupancyTestML")
195  << "[DTOccupancyTestML] ME: " << nEvtsName << " not found!" << endl;
196  }
197 
198  // Fill the global summary
199  // Check for entire sectors off and report them on the global summary
200  //FIXME: TODO
201 
202  if (writeRootFile)
203  ntuple->AutoSave("SaveSelf");
204 }
205 
207  LogVerbatim("DTDQM|DTMonitorClient|DTOccupancyTestML") << "[DTOccupancyTestML] endjob called!";
208  if (writeRootFile) {
209  rootFile->cd();
210  ntuple->Write();
211  rootFile->Close();
212  }
213 }
214 
215 // --------------------------------------------------
216 
217 void DTOccupancyTestML::bookHistos(DQMStore::IBooker& ibooker, const int wheelId, string folder, string histoTag) {
218  // Set the current folder
219  stringstream wheel;
220  wheel << wheelId;
221 
222  ibooker.setCurrentFolder(topFolder(true));
223 
224  // build the histo name
225  string histoName = histoTag + "_W" + wheel.str();
226 
227  LogVerbatim("DTDQM|DTMonitorClient|DTOccupancyTestML")
228  << "[DTOccupancyTestML]: booking wheel histo:" << histoName << " (tag " << histoTag
229  << ") in: " << topFolder(true) + "Wheel" + wheel.str() + "/" + folder << endl;
230 
231  string histoTitle = "Occupancy summary WHEEL: " + wheel.str();
232  if (tpMode) {
233  histoTitle = "TP Occupancy summary WHEEL: " + wheel.str();
234  }
235 
236  wheelHistos[wheelId] = ibooker.book2D(histoName, histoTitle, 12, 1, 13, 4, 1, 5);
237  wheelHistos[wheelId]->setBinLabel(1, "MB1", 2);
238  wheelHistos[wheelId]->setBinLabel(2, "MB2", 2);
239  wheelHistos[wheelId]->setBinLabel(3, "MB3", 2);
240  wheelHistos[wheelId]->setBinLabel(4, "MB4", 2);
241  wheelHistos[wheelId]->setAxisTitle("sector", 1);
242 }
243 
244 string DTOccupancyTestML::getMEName(string histoTag, const DTChamberId& chId) {
245  stringstream wheel;
246  wheel << chId.wheel();
247  stringstream station;
248  station << chId.station();
249  stringstream sector;
250  sector << chId.sector();
251 
252  string folderRoot =
253  topFolder(false) + "Wheel" + wheel.str() + "/Sector" + sector.str() + "/Station" + station.str() + "/";
254 
255  // build the histo name
256  string histoName = histoTag + "_W" + wheel.str() + "_St" + station.str() + "_Sec" + sector.str();
257 
258  string histoname = folderRoot + histoName;
259 
260  return histoname;
261 }
262 
263 int DTOccupancyTestML::getIntegral(TH2F* histo, int firstBinX, int lastBinX, int firstBinY, int lastBinY, bool doall) {
264  int sum = 0;
265  for (Int_t i = firstBinX; i < lastBinX + 1; i++) {
266  for (Int_t j = firstBinY; j < lastBinY + 1; j++) {
267  if (histo->GetBinContent(i, j) > 0) {
268  if (!doall)
269  return 1;
270  sum += histo->GetBinContent(i, j);
271  }
272  }
273  }
274 
275  return sum;
276 }
277 
279  const DTChamberId& chId,
280  float& chamberPercentage,
281  tensorflow::GraphDef* graphDef,
282  tensorflow::Session* session) {
283  LogTrace("DTDQM|DTMonitorClient|DTOccupancyTestML") << "--- Occupancy test for chamber: " << chId << endl;
284 
285  // Initialize counters
286  int totalLayers = 0;
287  int badLayers = 0;
288 
289  // Loop over the super layers
290  for (int superLayer = 1; superLayer <= 3; superLayer++) {
291  int binYlow = ((superLayer - 1) * 4) + 1;
292 
293  // Skip for non-existent super layers
294  if (chId.station() == 4 && superLayer == 2)
295  continue;
296 
297  // Loop over layers
298  for (int layer = 1; layer <= 4; layer++) {
299  DTLayerId layID(chId, superLayer, layer);
300  int firstWire = muonGeom->layer(layID)->specificTopology().firstChannel();
301  int nWires = muonGeom->layer(layID)->specificTopology().channels();
302  int binY = binYlow + (layer - 1);
303  std::vector<float> layerOccupancy(nWires);
304  int channelId = 0;
305 
306  // Loop over cells within a layer
307  for (int cell = firstWire; cell != (nWires + firstWire); cell++) {
308  double cellOccupancy = histo->GetBinContent(cell, binY);
309  layerOccupancy.at(channelId) = cellOccupancy;
310  channelId++;
311  }
312 
313  int targetSize = 47;
314  std::vector<float> reshapedLayerOccupancy = interpolateLayers(layerOccupancy, nWires, targetSize);
315 
316  // Scale occupancy
317  float maxOccupancyInLayer = *std::max_element(reshapedLayerOccupancy.begin(), reshapedLayerOccupancy.end());
318 
319  if (maxOccupancyInLayer != 0) {
320  for (auto& val : reshapedLayerOccupancy)
321  val /= maxOccupancyInLayer;
322  }
323 
324  // Define input
325  tensorflow::Tensor input(tensorflow::DT_FLOAT, {1, targetSize});
326  for (int i = 0; i < targetSize; i++)
327  input.matrix<float>()(0, i) = reshapedLayerOccupancy[i];
328 
329  std::vector<tensorflow::Tensor> outputs;
330  tensorflow::run(session, {{"input_cnn_input", input}}, {"output_cnn/Softmax"}, &outputs);
331 
332  totalLayers++;
333  bool isBad = outputs[0].matrix<float>()(0, 1) > 0.95;
334  if (isBad)
335  badLayers++;
336  }
337  }
338 
339  // Calculate a fraction of good layers
340  chamberPercentage = 1.0 - static_cast<float>(badLayers) / totalLayers;
341 
342  if (badLayers > 8)
343  return 3; // 3 SL
344  if (badLayers > 4)
345  return 2; // 2 SL
346  if (badLayers > 0)
347  return 1; // 1 SL
348  return 0;
349 }
350 
351 std::vector<float> DTOccupancyTestML::interpolateLayers(std::vector<float> const& inputs, int size, int targetSize) {
352  // Reshape layers with linear interpolation
353  int interpolationFloor = 0;
354  float interpolationPoint = 0.;
355  float step = static_cast<float>(size) / targetSize;
356  std::vector<float> reshapedInput(targetSize);
357 
358  int limitInLoop = inputs.size();
359  limitInLoop = std::min(limitInLoop, targetSize) - 1;
360  for (int i = 0; i < limitInLoop; i++) {
361  interpolationFloor = static_cast<int>(std::floor(interpolationPoint));
362  // Interpolating here
363  reshapedInput.at(i) =
364  (interpolationPoint - interpolationFloor) * (inputs[interpolationFloor + 1] - inputs[interpolationFloor]) +
365  inputs[interpolationFloor];
366  interpolationPoint = step + interpolationPoint;
367  }
368  return reshapedInput;
369 }
370 
371 string DTOccupancyTestML::topFolder(bool isBooking) const {
372  if (tpMode)
373  return string("DT/10-TestPulses/");
374  if (isBooking)
375  return string("DT/01-Digi/ML");
376  return string("DT/01-Digi/");
377 }
size
Write out results.
Log< level::Info, true > LogVerbatim
int station() const
Return the station number.
Definition: DTChamberId.h:42
void dqmEndJob(DQMStore::IBooker &, DQMStore::IGetter &) override
Endjob.
ESGetTokenH3DDVariant esConsumes(std::string const &Record, edm::ConsumesCollector &)
Definition: DeDxTools.cc:283
int getIntegral(TH2F *histo, int, int, int, int, bool)
std::string fullPath() const
Definition: FileInPath.cc:161
virtual void setCurrentFolder(std::string const &fullpath)
Definition: DQMStore.cc:36
virtual void setEntries(double nentries)
set # of entries
GraphDef * loadGraphDef(const std::string &pbFile)
Definition: TensorFlow.cc:120
void bookHistos(DQMStore::IBooker &, const int wheelId, std::string folder, std::string histoTag)
book the summary histograms
std::vector< float > interpolateLayers(std::vector< float > const &inputs, int size, int targetSize)
MonitorElement * summaryHisto
Definition: ntuple.py:1
#define LogTrace(id)
static std::string const input
Definition: EdmProvDump.cc:50
void dqmEndLuminosityBlock(DQMStore::IBooker &, DQMStore::IGetter &, edm::LuminosityBlock const &, edm::EventSetup const &) override
DQM Client Diagnostic.
T getUntrackedParameter(std::string const &, T const &) const
int firstChannel() const
Returns the wire number of the first wire.
Definition: DTTopology.h:79
void Fill(long long x)
virtual TH2F * getTH2F() const
virtual void Reset()
Remove all data from the ME, keept the empty histogram with all its settings.
std::string getMEName(std::string histoTag, const DTChamberId &chId)
Get the ME name.
const DTGeometry * muonGeom
virtual double getFloatValue() const
void run(Session *session, const NamedTensorList &inputs, const std::vector< std::string > &outputNames, std::vector< Tensor > *outputs, const thread::ThreadPoolOptions &threadPoolOptions)
Definition: TensorFlow.cc:259
bool closeSession(Session *&session)
Definition: TensorFlow.cc:234
void beginRun(edm::Run const &run, edm::EventSetup const &context) override
BeginRun.
Transition
Definition: Transition.h:12
MonitorElement * glbSummaryHisto
const DTTopology & specificTopology() const
Definition: DTLayer.cc:37
int runOccupancyTest(TH2F *histo, const DTChamberId &chId, float &chamberPercentage, tensorflow::GraphDef *graphDef, tensorflow::Session *session)
Session * createSession()
Definition: TensorFlow.cc:137
void setLogging(const std::string &level="3")
Definition: TensorFlow.cc:81
virtual void setBinContent(int binx, double content)
set content of bin (1-D)
MonitorElement * book2D(TString const &name, TString const &title, int nchX, double lowX, double highX, int nchY, double lowY, double highY, FUNC onbooking=NOOP())
Definition: DQMStore.h:212
virtual MonitorElement * get(std::string const &fullpath) const
Definition: DQMStore.cc:712
~DTOccupancyTestML() override
Destructor.
int wheel() const
Return the wheel number.
Definition: DTChamberId.h:39
HLT enums.
int sector() const
Definition: DTChamberId.h:49
edm::ESGetToken< DTGeometry, MuonGeometryRecord > muonGeomToken_
std::map< int, MonitorElement * > wheelHistos
std::string topFolder(bool isBooking) const
const std::vector< const DTChamber * > & chambers() const
Return a vector of all Chamber.
Definition: DTGeometry.cc:84
step
Definition: StallMonitor.cc:98
static char chambers[264][20]
Definition: ReadPGInfo.cc:243
std::string nameMonitoredHisto
DTOccupancyTestML(const edm::ParameterSet &ps)
Constructor.
int channels() const
Returns the number of wires in the layer.
Definition: DTTopology.h:76
Definition: Run.h:45
virtual double getBinContent(int binx) const
get content of bin (1-D)
const DTLayer * layer(const DTLayerId &id) const
Return a layer given its id.
Definition: DTGeometry.cc:96
virtual void setAxisTitle(const std::string &title, int axis=1)
set x-, y- or z-axis title (axis=1, 2, 3 respectively)