File indexing completed on 2021-11-12 00:00:59
0001 #include "DQM/EcalMonitorClient/interface/MLClient.h"
0002
0003 #include "DataFormats/EcalDetId/interface/EcalTrigTowerDetId.h"
0004
0005 #include "CondFormats/EcalObjects/interface/EcalDQMStatusHelper.h"
0006
0007 #include "DQM/EcalCommon/interface/EcalDQMCommonUtils.h"
0008
0009 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0010
0011 #include "PhysicsTools/ONNXRuntime/interface/ONNXRuntime.h"
0012
0013 #include "DQM/EcalCommon/interface/MESetNonObject.h"
0014
0015 using namespace cms::Ort;
0016
0017 namespace ecaldqm {
0018 MLClient::MLClient() : DQWorkerClient() { qualitySummaries_.insert("MLQualitySummary"); }
0019
0020 void MLClient::setParams(edm::ParameterSet const& _params) {
0021 MLThreshold_ = _params.getUntrackedParameter<double>("MLThreshold");
0022 PUcorr_slope_ = _params.getUntrackedParameter<double>("PUcorr_slope");
0023 PUcorr_intercept_ = _params.getUntrackedParameter<double>("PUcorr_intercept");
0024 avgOcc_ = _params.getUntrackedParameter<std::vector<double>>("avgOcc");
0025 if (!onlineMode_) {
0026 MEs_.erase(std::string("MLQualitySummary"));
0027 MEs_.erase(std::string("EventsperMLImage"));
0028 sources_.erase(std::string("PU"));
0029 sources_.erase(std::string("NumEvents"));
0030 sources_.erase(std::string("DigiAllByLumi"));
0031 sources_.erase(std::string("AELoss"));
0032 }
0033 }
0034
0035 void MLClient::producePlots(ProcessType) {
0036 if (!onlineMode_)
0037 return;
0038 using namespace std;
0039 MESet& meMLQualitySummary(MEs_.at("MLQualitySummary"));
0040 MESet& meEventsperMLImage(MEs_.at("EventsperMLImage"));
0041
0042 MESetNonObject const& sPU(static_cast<MESetNonObject&>(sources_.at("PU")));
0043 MESetNonObject const& sNumEvents(static_cast<MESetNonObject&>(sources_.at("NumEvents")));
0044
0045
0046 int nEv = sNumEvents.getFloatValue();
0047 double pu = sPU.getFloatValue();
0048
0049 if (pu < 0.) {
0050 return;
0051 }
0052 uint32_t mask(1 << EcalDQMStatusHelper::PEDESTAL_ONLINE_HIGH_GAIN_RMS_ERROR |
0053 1 << EcalDQMStatusHelper::PHYSICS_BAD_CHANNEL_WARNING |
0054 1 << EcalDQMStatusHelper::PHYSICS_BAD_CHANNEL_ERROR);
0055
0056
0057
0058
0059
0060
0061
0062
0063
0064
0065
0066
0067
0068
0069
0070
0071 TH2F* hEbDigiMap((sources_.at("DigiAllByLumi")).getME(1)->getTH2F());
0072
0073 size_t nTowers = nEtaTowers * nPhiTowers;
0074 std::vector<float> ebOccMap1dCumulPad;
0075 std::valarray<float> ebOccMap1d(nTowers);
0076
0077
0078 for (int i = 0; i < hEbDigiMap->GetNbinsY(); i++) {
0079 for (int j = 0; j < hEbDigiMap->GetNbinsX(); j++) {
0080 int bin = hEbDigiMap->GetBin(j + 1, i + 1);
0081 int k = (i * nPhiTowers) + j;
0082 ebOccMap1d[k] = hEbDigiMap->GetBinContent(bin);
0083 }
0084 }
0085 ebOccMap1dQ.push_back(ebOccMap1d);
0086 NEventQ.push_back(nEv);
0087
0088 if (NEventQ.size() < nLS) {
0089 return;
0090 }
0091 if (NEventQ.size() > nLS) {
0092 NEventQ.pop_front();
0093 }
0094 if (ebOccMap1dQ.size() > nLS) {
0095 ebOccMap1dQ.pop_front();
0096 }
0097
0098 int TNum = 0;
0099 for (size_t i = 0; i < nLS; i++) {
0100 TNum += NEventQ[i];
0101 }
0102 if (TNum < 200) {
0103 return;
0104 }
0105
0106 meEventsperMLImage.fill(getEcalDQMSetupObjects(), EcalBarrel, double(timestamp_.iLumi), double(TNum));
0107
0108
0109 std::valarray<float> ebOccMap1dCumul(0., nTowers);
0110
0111 for (size_t i = 0; i < ebOccMap1dQ.size(); i++) {
0112 ebOccMap1dCumul += ebOccMap1dQ[i];
0113 }
0114
0115 ebOccMap1dCumul = ebOccMap1dCumul / (PUcorr_slope_ * pu + PUcorr_intercept_);
0116
0117
0118
0119 ebOccMap1dCumul = ebOccMap1dCumul * (nEtaTowersPad * nPhiTowers);
0120
0121
0122 ebOccMap1dCumul = ebOccMap1dCumul * (500. / TNum);
0123
0124
0125 ebOccMap1dCumulPad.assign(std::begin(ebOccMap1dCumul), std::end(ebOccMap1dCumul));
0126
0127
0128 for (int k = 0; k < nPhiTowers; k++) {
0129 float val = ebOccMap1dCumulPad[nPhiTowers - 1];
0130 ebOccMap1dCumulPad.insert(ebOccMap1dCumulPad.begin(),
0131 val);
0132 }
0133
0134 int size = ebOccMap1dCumulPad.size();
0135 for (int k = (size - nPhiTowers); k < size; k++) {
0136 float val = ebOccMap1dCumulPad[k];
0137 ebOccMap1dCumulPad.push_back(val);
0138 }
0139
0140
0141
0142
0143
0144
0145
0146
0147
0148
0149
0150
0151
0152
0153
0154
0155
0156
0157 std::string instanceName{"AE-DQM-inference"};
0158 std::string modelFilepath = edm::FileInPath("DQM/EcalMonitorClient/data/onnxModels/resnet.onnx").fullPath();
0159
0160 Ort::SessionOptions sessionOptions;
0161 sessionOptions.SetIntraOpNumThreads(1);
0162 Ort::Env env(OrtLoggingLevel::ORT_LOGGING_LEVEL_WARNING, instanceName.c_str());
0163 Ort::Session session(env, modelFilepath.c_str(), sessionOptions);
0164
0165 Ort::AllocatorWithDefaultOptions allocator;
0166
0167 const char* inputName = session.GetInputName(0, allocator);
0168
0169 Ort::TypeInfo inputTypeInfo = session.GetInputTypeInfo(0);
0170 auto inputTensorInfo = inputTypeInfo.GetTensorTypeAndShapeInfo();
0171
0172 std::vector<int64_t> inputDims = inputTensorInfo.GetShape();
0173
0174 const char* outputName = session.GetOutputName(0, allocator);
0175
0176 Ort::TypeInfo outputTypeInfo = session.GetOutputTypeInfo(0);
0177 auto outputTensorInfo = outputTypeInfo.GetTensorTypeAndShapeInfo();
0178
0179 std::vector<int64_t> outputDims = outputTensorInfo.GetShape();
0180
0181 size_t TensorSize = nEtaTowersPad * nPhiTowers;
0182 std::vector<float> ebRecoOccMap1dPad(TensorSize);
0183
0184 std::vector<const char*> inputNames{inputName};
0185 std::vector<const char*> outputNames{outputName};
0186 std::vector<Ort::Value> inputTensors;
0187 std::vector<Ort::Value> outputTensors;
0188
0189 Ort::MemoryInfo memoryInfo =
0190 Ort::MemoryInfo::CreateCpu(OrtAllocatorType::OrtArenaAllocator, OrtMemType::OrtMemTypeDefault);
0191 inputTensors.push_back(Ort::Value::CreateTensor<float>(
0192 memoryInfo, ebOccMap1dCumulPad.data(), TensorSize, inputDims.data(), inputDims.size()));
0193
0194 outputTensors.push_back(Ort::Value::CreateTensor<float>(
0195 memoryInfo, ebRecoOccMap1dPad.data(), TensorSize, outputDims.data(), outputDims.size()));
0196
0197 session.Run(Ort::RunOptions{nullptr},
0198 inputNames.data(),
0199 inputTensors.data(),
0200 1,
0201 outputNames.data(),
0202 outputTensors.data(),
0203 1);
0204
0205
0206
0207
0208 std::valarray<std::valarray<float>> lossMap2d(std::valarray<float>(nPhiTowers), nEtaTowers);
0209
0210
0211
0212 std::valarray<float> recoOcc1d(0., nPhiTowers);
0213 std::valarray<float> inputOcc1d(0., nPhiTowers);
0214 std::valarray<float> avgOcc1d(0., nPhiTowers);
0215 std::valarray<float> loss_;
0216
0217
0218
0219
0220 for (int i = 1; i < 35; i++) {
0221 for (int j = 0; j < nPhiTowers; j++) {
0222 int k = (i * nPhiTowers) + j;
0223 recoOcc1d[j] = ebRecoOccMap1dPad[k];
0224 inputOcc1d[j] = ebOccMap1dCumulPad[k];
0225 avgOcc1d[j] = avgOcc_[k];
0226 }
0227
0228 loss_ = std::pow((recoOcc1d / avgOcc1d - inputOcc1d / avgOcc1d), 2);
0229 lossMap2d[i - 1] = (loss_);
0230 }
0231
0232 lossMap2dQ.push_back(lossMap2d);
0233 if (lossMap2dQ.size() > nLSloss) {
0234 lossMap2dQ.pop_front();
0235 }
0236 if (lossMap2dQ.size() < nLSloss) {
0237 return;
0238 }
0239
0240 std::valarray<std::valarray<float>> lossMap2dMult(std::valarray<float>(1., nPhiTowers), nEtaTowers);
0241
0242
0243
0244 for (size_t i = 0; i < lossMap2dQ.size(); i++) {
0245 lossMap2dMult *= lossMap2dQ[i];
0246 }
0247
0248
0249 MESet const& sAELoss(sources_.at("AELoss"));
0250 TH2F* hLossMap2dMult(sAELoss.getME(1)->getTH2F());
0251 for (int i = 0; i < hLossMap2dMult->GetNbinsY(); i++) {
0252 for (int j = 0; j < hLossMap2dMult->GetNbinsX(); j++) {
0253 int bin_ = hLossMap2dMult->GetBin(j + 1, i + 1);
0254 double content = lossMap2dMult[i][j];
0255 hLossMap2dMult->SetBinContent(bin_, content);
0256 }
0257 }
0258
0259
0260
0261
0262 MESet::const_iterator dAEnd(sAELoss.end(GetElectronicsMap()));
0263 for (MESet::const_iterator dItr(sAELoss.beginChannel(GetElectronicsMap())); dItr != dAEnd;
0264 dItr.toNextChannel(GetElectronicsMap())) {
0265 DetId id(dItr->getId());
0266
0267 bool doMaskML(meMLQualitySummary.maskMatches(id, mask, statusManager_, GetTrigTowerMap()));
0268
0269 float entries(dItr->getBinContent());
0270 int quality(doMaskML ? kMGood : kGood);
0271
0272 if (entries > MLThreshold_) {
0273 quality = doMaskML ? kMBad : kBad;
0274 }
0275
0276 meMLQualitySummary.setBinContent(getEcalDQMSetupObjects(), id, double(quality));
0277 }
0278 }
0279
0280 DEFINE_ECALDQM_WORKER(MLClient);
0281 }