EcalMonitorClient/src/MLClient.cc

0001 #include "DQM/EcalMonitorClient/interface/MLClient.h"
0002
0003 #include "CondFormats/EcalObjects/interface/EcalDQMStatusHelper.h"
0004
0005 #include "DQM/EcalCommon/interface/EcalDQMCommonUtils.h"
0006
0007 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0008
0009 #include "PhysicsTools/ONNXRuntime/interface/ONNXRuntime.h"
0010
0011 #include "DQM/EcalCommon/interface/MESetNonObject.h"
0012
0013 #include <fstream>
0014 #include <string>
0015
0016 using namespace cms::Ort;
0017
0018 namespace ecaldqm {
0019   MLClient::MLClient() : DQWorkerClient() { qualitySummaries_.insert("MLQualitySummary"); }
0020
0021   void MLClient::setParams(edm::ParameterSet const& _params) {
0022     EBThreshold_ = _params.getUntrackedParameter<double>("EBThreshold");
0023     EEpThreshold_ = _params.getUntrackedParameter<double>("EEpThreshold");
0024     EEmThreshold_ = _params.getUntrackedParameter<double>("EEmThreshold");
0025     EB_PUcorr_slope_ = _params.getUntrackedParameter<double>("EB_PUcorr_slope");
0026     EB_PUcorr_intercept_ = _params.getUntrackedParameter<double>("EB_PUcorr_intercept");
0027     EEp_PUcorr_slope_ = _params.getUntrackedParameter<double>("EEp_PUcorr_slope");
0028     EEp_PUcorr_intercept_ = _params.getUntrackedParameter<double>("EEp_PUcorr_intercept");
0029     EEm_PUcorr_slope_ = _params.getUntrackedParameter<double>("EEm_PUcorr_slope");
0030     EEm_PUcorr_intercept_ = _params.getUntrackedParameter<double>("EEm_PUcorr_intercept");
0031
0032     if (!onlineMode_) {
0033       MEs_.erase(std::string("MLQualitySummary"));
0034       MEs_.erase(std::string("EventsperMLImage"));
0035       sources_.erase(std::string("PU"));
0036       sources_.erase(std::string("NumEvents"));
0037       sources_.erase(std::string("DigiAllByLumi"));
0038       sources_.erase(std::string("AELoss"));
0039       sources_.erase(std::string("BadTowerCount"));
0040       sources_.erase(std::string("BadTowerCountNorm"));
0041     }
0042   }
0043
0044   void MLClient::producePlots(ProcessType) {
0045     if (!onlineMode_)
0046       return;
0047     nbadtowerEB = 0;
0048     nbadtowerEE = 0;
0049
0050     using namespace std;
0051     MESet& meMLQualitySummary(MEs_.at("MLQualitySummary"));
0052     MESet& meEventsperMLImage(MEs_.at("EventsperMLImage"));
0053
0054     MESetNonObject const& sPU(static_cast<MESetNonObject&>(sources_.at("PU")));
0055     MESetNonObject const& sNumEvents(static_cast<MESetNonObject&>(sources_.at("NumEvents")));
0056
0057     //Get the no.of events and the PU per LS calculated in OccupancyTask
0058     int nEv = sNumEvents.getFloatValue();
0059     double pu = sPU.getFloatValue();
0060
0061     //Do not compute ML quality if PU is non existent.
0062     if (pu <= 0.) {
0063       return;
0064     }
0065     uint32_t mask(1 << EcalDQMStatusHelper::PEDESTAL_ONLINE_HIGH_GAIN_RMS_ERROR |
0066                   1 << EcalDQMStatusHelper::PHYSICS_BAD_CHANNEL_WARNING |
0067                   1 << EcalDQMStatusHelper::PHYSICS_BAD_CHANNEL_ERROR);
0068
0069     //////////////// ML Data Preprocessing //////////////////////////////////
0070     //Inorder to feed the data into the ML model we apply some preprocessing.
0071     //We use the Digi Occupancy per Lumisection as the input source.
0072     //The model was trained on each occupancy plot having 500 events.
0073     //In apprehension of the low luminosity in the beginning of Run3, where in online DQM
0074     //the no.of events per LS could be lower than 500, we sum the occupancies over a fixed no.of lumisections as a running sum,
0075     //and require that the total no.of events on this summed occupancy to be atleast 200.
0076     //(This no.of LS and the no.of events are parameters which would require tuning later)
0077     //This summed occupancy is now the input image, which is then corrected for PileUp(PU) dependence and
0078     //change in no.of events, which are derived from training.
0079     //The input image is also padded by replicating the top and bottom rows so as to prevent the "edge effect"
0080     //wherein the ML model's learning degrades near the edge of the data set it sees.
0081     //This padding is then removed during inference on the model output.
0082
0083     //Get the histogram of the input digi occupancy per lumisection.
0084     TH2F* hEEmDigiMap((sources_.at("DigiAllByLumi")).getME(0)->getTH2F());
0085     TH2F* hEbDigiMap((sources_.at("DigiAllByLumi")).getME(1)->getTH2F());
0086     TH2F* hEEpDigiMap((sources_.at("DigiAllByLumi")).getME(2)->getTH2F());
0087
0088     size_t nEBTowers = nEBEtaTowers * nEBPhiTowers;  //Each EB occupancy map is of size 34x72 towers
0089     size_t nEETowers = nEEEtaTowers * nEEPhiTowers;  //Each EE occupancy map is of size 20x20 towers
0090
0091     //Vectors to feed into the ML network
0092     std::vector<float> ebOccMap1dCumulPad;
0093     std::vector<float> eemOccMap1dCumulPad;
0094     std::vector<float> eepOccMap1dCumulPad;
0095
0096     //Array to store occupancy maps
0097     std::valarray<float> ebOccMap1d(nEBTowers);
0098     std::valarray<float> eemOccMap1d(nEETowers);
0099     std::valarray<float> eepOccMap1d(nEETowers);
0100
0101     //Store the values from the input histogram into the array
0102     //to do preprocessing
0103     for (int i = 0; i < hEbDigiMap->GetNbinsY(); i++) {  //NbinsY = 34, NbinsX = 72
0104       for (int j = 0; j < hEbDigiMap->GetNbinsX(); j++) {
0105         int bin = hEbDigiMap->GetBin(j + 1, i + 1);
0106         int k = (i * nEBPhiTowers) + j;
0107         ebOccMap1d[k] = hEbDigiMap->GetBinContent(bin);
0108       }
0109     }
0110     ebOccMap1dQ.push_back(ebOccMap1d);  //Queue which stores input occupancy maps for nLS lumis
0111
0112     for (int i = 0; i < hEEpDigiMap->GetNbinsY(); i++) {  //NbinsY = 20, NbinsX = 20
0113       for (int j = 0; j < hEEpDigiMap->GetNbinsX(); j++) {
0114         int bin = hEEpDigiMap->GetBin(j + 1, i + 1);
0115         int k = (i * nEEPhiTowers) + j;
0116         eemOccMap1d[k] = hEEmDigiMap->GetBinContent(bin);
0117         eepOccMap1d[k] = hEEpDigiMap->GetBinContent(bin);
0118       }
0119     }
0120
0121     //Queue which stores input occupancy maps for nLS lumis
0122     eemOccMap1dQ.push_back(eemOccMap1d);
0123     eepOccMap1dQ.push_back(eepOccMap1d);
0124
0125     NEventQ.push_back(nEv);  //Queue which stores the no.of events per LS for nLS lumis
0126
0127     if (NEventQ.size() < nLS) {
0128       return;  //Should have nLS lumis to add the occupancy over.
0129     }
0130     if (NEventQ.size() > nLS) {
0131       NEventQ.pop_front();  //Keep only nLS consecutive LS. Pop the first one if size greater than nLS
0132     }
0133     if (ebOccMap1dQ.size() > nLS) {
0134       ebOccMap1dQ.pop_front();  //Same conditon for the input occupancy maps.
0135       eemOccMap1dQ.pop_front();
0136       eepOccMap1dQ.pop_front();
0137     }
0138
0139     int TNum = 0;
0140     for (size_t i = 0; i < nLS; i++) {
0141       TNum += NEventQ[i];  //Total no.of events over nLS lumis
0142     }
0143
0144     if (TNum < 400) {
0145       return;  //The total no.of events should be atleast 400 over nLS for meaningful statistics
0146     }
0147     //Fill the ME to monitor the trend of the total no.of events in each input image to the ML model
0148     meEventsperMLImage.fill(getEcalDQMSetupObjects(), EcalBarrel, double(timestamp_.iLumi), double(TNum));
0149
0150     //Array to hold the sum of inputs, which make atleast 400 events.
0151     std::valarray<float> ebOccMap1dCumul(0., nEBTowers);
0152     std::valarray<float> eemOccMap1dCumul(0., nEETowers);
0153     std::valarray<float> eepOccMap1dCumul(0., nEETowers);
0154
0155     //Sum the input arrays of nLS.
0156     for (size_t i = 0; i < ebOccMap1dQ.size(); i++) {
0157       ebOccMap1dCumul += ebOccMap1dQ[i];
0158       eemOccMap1dCumul += eemOccMap1dQ[i];
0159       eepOccMap1dCumul += eepOccMap1dQ[i];
0160     }
0161
0162     //Applying PU correction derived from training
0163     ebOccMap1dCumul = ebOccMap1dCumul / (EB_PUcorr_slope_ * pu + EB_PUcorr_intercept_);
0164     eemOccMap1dCumul = eemOccMap1dCumul / (EEm_PUcorr_slope_ * pu + EEm_PUcorr_intercept_);
0165     eepOccMap1dCumul = eepOccMap1dCumul / (EEp_PUcorr_slope_ * pu + EEp_PUcorr_intercept_);
0166
0167     //Scaling up to match input dimensions.
0168     ebOccMap1dCumul = ebOccMap1dCumul * (nEBEtaTowers * nEBPhiTowers);
0169     eemOccMap1dCumul = eemOccMap1dCumul * nEEEtaTowers * nEEPhiTowers;  //(nEETowersPad * nEETowersPad);
0170     eepOccMap1dCumul = eepOccMap1dCumul * nEEEtaTowers * nEEPhiTowers;  //(nEETowersPad * nEETowersPad);
0171
0172     //Correction for no.of events in each input image as originally model trained with 500 events per image
0173     ebOccMap1dCumul = ebOccMap1dCumul * (500. / TNum);
0174     eemOccMap1dCumul = eemOccMap1dCumul * (500. / TNum);
0175     eepOccMap1dCumul = eepOccMap1dCumul * (500. / TNum);
0176
0177     std::vector<std::vector<float>> ebOccMap2dCumul(nEBEtaTowers, std::vector<float>(nEBPhiTowers, 0.));
0178     //Convert 1dCumul to 2d
0179     for (size_t i = 0; i < nEBEtaTowers; i++) {
0180       for (size_t j = 0; j < nEBPhiTowers; j++) {
0181         int k = (i * nEBPhiTowers) + j;
0182         ebOccMap2dCumul[i][j] = ebOccMap1dCumul[k];
0183       }
0184     }
0185
0186     std::vector<float> pad_top;
0187     std::vector<float> pad_bottom;
0188     std::vector<float> pad_left;
0189     std::vector<float> pad_right;
0190
0191     pad_top = ebOccMap2dCumul[0];
0192     pad_bottom = ebOccMap2dCumul[ebOccMap2dCumul.size() - 1];
0193
0194     ebOccMap2dCumul.insert(ebOccMap2dCumul.begin(), pad_top);
0195     ebOccMap2dCumul.push_back(pad_bottom);
0196
0197     //// Endcaps ///
0198     std::vector<std::vector<float>> eemOccMap2dCumul(nEEEtaTowers, std::vector<float>(nEEPhiTowers, 0.));
0199     std::vector<std::vector<float>> eepOccMap2dCumul(nEEEtaTowers, std::vector<float>(nEEPhiTowers, 0.));
0200
0201     for (size_t i = 0; i < nEEEtaTowers; i++) {
0202       for (size_t j = 0; j < nEEPhiTowers; j++) {
0203         int k = (i * nEEPhiTowers) + j;
0204         eemOccMap2dCumul[i][j] = eemOccMap1dCumul[k];
0205         eepOccMap2dCumul[i][j] = eepOccMap1dCumul[k];
0206       }
0207     }
0208
0209     // EE - //
0210     pad_top.clear();
0211     pad_bottom.clear();
0212     pad_left.clear();
0213     pad_right.clear();
0214
0215     pad_top = eemOccMap2dCumul[0];
0216     pad_bottom = eemOccMap2dCumul[eemOccMap2dCumul.size() - 1];
0217
0218     eemOccMap2dCumul.insert(eemOccMap2dCumul.begin(), pad_top);
0219     eemOccMap2dCumul.push_back(pad_bottom);
0220
0221     for (auto& row : eemOccMap2dCumul) {
0222       pad_left.push_back(row[0]);
0223       pad_right.push_back(row[row.size() - 1]);
0224     }
0225
0226     std::size_t Lindex = 0;
0227     std::size_t Rindex = 0;
0228
0229     for (auto& row : eemOccMap2dCumul) {
0230       row.insert(row.begin(), pad_left[Lindex++]);
0231       row.insert(row.end(), pad_right[Rindex++]);
0232     }
0233
0234     // EE + //
0235     pad_top.clear();
0236     pad_bottom.clear();
0237
0238     pad_top = eepOccMap2dCumul[0];
0239     pad_bottom = eepOccMap2dCumul[eepOccMap2dCumul.size() - 1];
0240
0241     eepOccMap2dCumul.insert(eepOccMap2dCumul.begin(), pad_top);
0242     eepOccMap2dCumul.push_back(pad_bottom);
0243
0244     for (auto& row : eepOccMap2dCumul) {
0245       pad_left.push_back(row[0]);
0246       pad_right.push_back(row[row.size() - 1]);
0247     }
0248
0249     Lindex = 0;
0250     Rindex = 0;
0251
0252     for (auto& row : eepOccMap2dCumul) {
0253       row.insert(row.begin(), pad_left[Lindex++]);
0254       row.insert(row.end(), pad_right[Rindex++]);
0255     }
0256
0257     //The pre-processed input is now fed into the 1D input tensor vector which will go into the ML model
0258     for (auto& row : ebOccMap2dCumul) {
0259       ebOccMap1dCumulPad.insert(ebOccMap1dCumulPad.end(), row.begin(), row.end());
0260     }
0261
0262     for (auto& row : eemOccMap2dCumul) {
0263       eemOccMap1dCumulPad.insert(eemOccMap1dCumulPad.end(), row.begin(), row.end());
0264     }
0265
0266     for (auto& row : eepOccMap2dCumul) {
0267       eepOccMap1dCumulPad.insert(eepOccMap1dCumulPad.end(), row.begin(), row.end());
0268     }
0269
0270     ///// Model Inference //////
0271     //An Autoencoder (AE) network with resnet architecture is used here which is trained on
0272     //certified good data (EB digi occupancy) from Run 2018 data.
0273     //On giving an input occupancy map, the encoder part of the AE compresses and reduces the input data, learning its features,
0274     //and the decoder reconstructs the data from the encoded form into a representation as close to the original input as possible.
0275     //We then compute the Mean squared error (MSE) between the input and output image, also called the Reconstruction Loss,
0276     //calculated at a tower by tower basis.
0277     //Thus, given an anomalous tower the loss should be significantly higher than the loss with respect to good towers, which the model
0278     //has already seen --> anomaly detection.
0279     //When calculating the loss we also apply a response correction by dividing each input and output image with the average occupancy from
0280     //all 2018 data (also to be tuned),to accommodate the difference in response of crystals in different regions of the Ecal Barrel
0281     //Further each loss map from each input image is then multiplied by the last N loss maps,
0282     ///so that real anomalies which persist with time are enhanced and fluctuations are suppressed.
0283     //A quality threshold is then applied on this time multiplied loss map, to mark them as GOOD or BAD,
0284     //after which it is stored as a quality summary ME.
0285
0286     ///ONNX model running///
0287     std::string instanceName{"AE-DQM-inference"};
0288     std::string modelFilepath = edm::FileInPath("DQM/EcalMonitorClient/data/onnxModels/resnet.onnx").fullPath();
0289
0290     Ort::SessionOptions sessionOptions;
0291     sessionOptions.SetIntraOpNumThreads(1);
0292     Ort::Env env(OrtLoggingLevel::ORT_LOGGING_LEVEL_WARNING, instanceName.c_str());
0293     Ort::Session session(env, modelFilepath.c_str(), sessionOptions);
0294
0295     Ort::AllocatorWithDefaultOptions allocator;
0296
0297     // Strings returned by session.GetInputNameAllocated are temporary, need to copy them before they are deallocated
0298     std::string inputName{session.GetInputNameAllocated(0, allocator).get()};
0299
0300     Ort::TypeInfo inputTypeInfo = session.GetInputTypeInfo(0);
0301     auto inputTensorInfo = inputTypeInfo.GetTensorTypeAndShapeInfo();
0302
0303     std::vector<int64_t> inputDims = inputTensorInfo.GetShape();
0304
0305     std::string outputName{session.GetOutputNameAllocated(0, allocator).get()};
0306
0307     Ort::TypeInfo outputTypeInfo = session.GetOutputTypeInfo(0);
0308     auto outputTensorInfo = outputTypeInfo.GetTensorTypeAndShapeInfo();
0309
0310     std::vector<int64_t> outputDims = outputTensorInfo.GetShape();
0311
0312     size_t TensorSize = nEBEtaTowersPad * nEBPhiTowers;
0313     std::vector<float> ebRecoOccMap1dPad(TensorSize);  //To store the output reconstructed occupancy
0314
0315     std::vector<const char*> inputNames{inputName.c_str()};
0316     std::vector<const char*> outputNames{outputName.c_str()};
0317     std::vector<Ort::Value> inputTensors;
0318     std::vector<Ort::Value> outputTensors;
0319
0320     Ort::MemoryInfo memoryInfo =
0321         Ort::MemoryInfo::CreateCpu(OrtAllocatorType::OrtArenaAllocator, OrtMemType::OrtMemTypeDefault);
0322     inputTensors.push_back(Ort::Value::CreateTensor<float>(
0323         memoryInfo, ebOccMap1dCumulPad.data(), TensorSize, inputDims.data(), inputDims.size()));
0324
0325     outputTensors.push_back(Ort::Value::CreateTensor<float>(
0326         memoryInfo, ebRecoOccMap1dPad.data(), TensorSize, outputDims.data(), outputDims.size()));
0327
0328     session.Run(Ort::RunOptions{nullptr},
0329                 inputNames.data(),
0330                 inputTensors.data(),
0331                 1,
0332                 outputNames.data(),
0333                 outputTensors.data(),
0334                 1);
0335
0336     //Endcaps
0337     // EE- //
0338
0339     inputDims.clear();
0340     outputDims.clear();
0341     inputNames.clear();
0342     outputNames.clear();
0343     inputTensors.clear();
0344     outputTensors.clear();
0345
0346     modelFilepath = edm::FileInPath("DQM/EcalMonitorClient/data/onnxModels/EEm_resnet2018.onnx").fullPath();
0347
0348     Ort::Session EEm_session(env, modelFilepath.c_str(), sessionOptions);
0349
0350     inputName = EEm_session.GetInputNameAllocated(0, allocator).get();
0351
0352     inputTypeInfo = EEm_session.GetInputTypeInfo(0);
0353     auto EEm_inputTensorInfo = inputTypeInfo.GetTensorTypeAndShapeInfo();
0354
0355     inputDims = EEm_inputTensorInfo.GetShape();
0356
0357     outputName = EEm_session.GetOutputNameAllocated(0, allocator).get();
0358
0359     //Ort::TypeInfo
0360     outputTypeInfo = EEm_session.GetOutputTypeInfo(0);
0361     auto EEm_outputTensorInfo = outputTypeInfo.GetTensorTypeAndShapeInfo();
0362
0363     outputDims = EEm_outputTensorInfo.GetShape();
0364
0365     size_t EE_TensorSize = nEETowersPad * nEETowersPad;
0366     std::vector<float> eemRecoOccMap1dPad(EE_TensorSize);  //To store the output reconstructed occupancy
0367
0368     inputNames.push_back(inputName.c_str());
0369     outputNames.push_back(outputName.c_str());
0370
0371     //Ort::MemoryInfo
0372     memoryInfo = Ort::MemoryInfo::CreateCpu(OrtAllocatorType::OrtArenaAllocator, OrtMemType::OrtMemTypeDefault);
0373     inputTensors.push_back(Ort::Value::CreateTensor<float>(
0374         memoryInfo, eemOccMap1dCumulPad.data(), EE_TensorSize, inputDims.data(), inputDims.size()));
0375
0376     outputTensors.push_back(Ort::Value::CreateTensor<float>(
0377         memoryInfo, eemRecoOccMap1dPad.data(), EE_TensorSize, outputDims.data(), outputDims.size()));
0378
0379     EEm_session.Run(Ort::RunOptions{nullptr},
0380                     inputNames.data(),
0381                     inputTensors.data(),
0382                     1,
0383                     outputNames.data(),
0384                     outputTensors.data(),
0385                     1);
0386
0387     // EE+ //
0388     inputDims.clear();
0389     outputDims.clear();
0390     inputNames.clear();
0391     outputNames.clear();
0392     inputTensors.clear();
0393     outputTensors.clear();
0394
0395     modelFilepath = edm::FileInPath("DQM/EcalMonitorClient/data/onnxModels/EEp_resnet2018.onnx").fullPath();
0396
0397     Ort::Session EEp_session(env, modelFilepath.c_str(), sessionOptions);
0398
0399     inputName = EEp_session.GetInputNameAllocated(0, allocator).get();
0400
0401     inputTypeInfo = EEp_session.GetInputTypeInfo(0);
0402     auto EEp_inputTensorInfo = inputTypeInfo.GetTensorTypeAndShapeInfo();
0403
0404     inputDims = EEp_inputTensorInfo.GetShape();
0405
0406     outputName = EEp_session.GetOutputNameAllocated(0, allocator).get();
0407
0408     outputTypeInfo = EEp_session.GetOutputTypeInfo(0);
0409
0410     auto EEp_outputTensorInfo = outputTypeInfo.GetTensorTypeAndShapeInfo();
0411
0412     outputDims = EEp_outputTensorInfo.GetShape();
0413
0414     std::vector<float> eepRecoOccMap1dPad(EE_TensorSize);  //To store the output reconstructed occupancy
0415
0416     inputNames.push_back(inputName.c_str());
0417     outputNames.push_back(outputName.c_str());
0418
0419     //Ort::MemoryInfo
0420     memoryInfo = Ort::MemoryInfo::CreateCpu(OrtAllocatorType::OrtArenaAllocator, OrtMemType::OrtMemTypeDefault);
0421     inputTensors.push_back(Ort::Value::CreateTensor<float>(
0422         memoryInfo, eepOccMap1dCumulPad.data(), EE_TensorSize, inputDims.data(), inputDims.size()));
0423
0424     outputTensors.push_back(Ort::Value::CreateTensor<float>(
0425         memoryInfo, eepRecoOccMap1dPad.data(), EE_TensorSize, outputDims.data(), outputDims.size()));
0426
0427     EEp_session.Run(Ort::RunOptions{nullptr},
0428                     inputNames.data(),
0429                     inputTensors.data(),
0430                     1,
0431                     outputNames.data(),
0432                     outputTensors.data(),
0433                     1);
0434
0435     ///Inference on the output from the model///
0436     //2D Loss map to store tower by tower loss between the output (reconstructed) and input occupancies,
0437     //Have same dimensions as the occupancy plot
0438     std::valarray<std::valarray<float>> EBlossMap2d(std::valarray<float>(nEBPhiTowers), nEBEtaTowers);
0439     std::valarray<std::valarray<float>> EEmlossMap2d(std::valarray<float>(nEEPhiTowers), nEEEtaTowers);
0440     std::valarray<std::valarray<float>> EEplossMap2d(std::valarray<float>(nEEPhiTowers), nEEEtaTowers);
0441
0442     //1D val arrays to store row wise information corresponding to the reconstructed, input and average occupancies, and loss.
0443     //and to do element wise (tower wise) operations on them to calculate the MSE loss between the reco and input occupancy.
0444     std::valarray<float> EBrecoOcc1d(0., nEBPhiTowers);
0445     std::valarray<float> EBinputOcc1d(0., nEBPhiTowers);
0446     std::valarray<float> EBavgOcc1d(0., nEBPhiTowers);
0447     std::valarray<float> EBloss_;
0448
0449     std::valarray<float> EEmrecoOcc1d(0., nEEPhiTowers);
0450     std::valarray<float> EEminputOcc1d(0., nEEPhiTowers);
0451     std::valarray<float> EEmavgOcc1d(0., nEEPhiTowers);
0452     std::valarray<float> EEmloss_;
0453
0454     std::valarray<float> EEprecoOcc1d(0., nEEPhiTowers);
0455     std::valarray<float> EEpinputOcc1d(0., nEEPhiTowers);
0456     std::valarray<float> EEpavgOcc1d(0., nEEPhiTowers);
0457     std::valarray<float> EEploss_;
0458
0459     std::string EBOccpath =
0460         edm::FileInPath("DQM/EcalMonitorClient/data/MLAvgOccupancy/EB_avgocc_Run2022_500ev.dat").fullPath();
0461     std::ifstream inFile;
0462     double val;
0463     inFile.open((EBOccpath).c_str());
0464     while (inFile) {
0465       inFile >> val;
0466       if (inFile.eof())
0467         break;
0468       EBavgOcc.push_back(val);
0469     }
0470     inFile.close();
0471
0472     std::string EEmOccpath =
0473         edm::FileInPath("DQM/EcalMonitorClient/data/MLAvgOccupancy/EEm_avgocc_Run2022_500ev.dat").fullPath();
0474     inFile.open((EEmOccpath).c_str());
0475     while (inFile) {
0476       inFile >> val;
0477       if (inFile.eof())
0478         break;
0479       EEmavgOcc.push_back(val);
0480     }
0481     inFile.close();
0482
0483     std::string EEpOccpath =
0484         edm::FileInPath("DQM/EcalMonitorClient/data/MLAvgOccupancy/EEp_avgocc_Run2022_500ev.dat").fullPath();
0485     inFile.open((EEpOccpath).c_str());
0486     while (inFile) {
0487       inFile >> val;
0488       if (inFile.eof())
0489         break;
0490       EEpavgOcc.push_back(val);
0491     }
0492     inFile.close();
0493
0494     //Loss calculation
0495     //Ignore the top and bottom replicated padded rows when doing inference
0496     //by making index i run over (1,35) instead of (0,36) for EB, and over (1,21) for EE
0497
0498     MESet const& sAEReco(sources_.at("AEReco"));
0499     TH2F* hEBRecoMap2d(sAEReco.getME(1)->getTH2F());
0500
0501     for (int i = 1; i < nEBEtaTowersPad - 1; i++) {
0502       for (int j = 0; j < nEBPhiTowers; j++) {
0503         int k = (i * nEBPhiTowers) + j;
0504         int bin_ = hEBRecoMap2d->GetBin(j + 1, i);
0505         EBrecoOcc1d[j] = ebRecoOccMap1dPad[k];
0506         EBinputOcc1d[j] = ebOccMap1dCumulPad[k];
0507         EBavgOcc1d[j] = EBavgOcc[k];
0508         double content = ebRecoOccMap1dPad[k];
0509         hEBRecoMap2d->SetBinContent(bin_, content);
0510       }
0511       //Calculate the MSE loss = (output-input)^2, with avg response correction
0512       EBloss_ = std::pow((EBrecoOcc1d / EBavgOcc1d - EBinputOcc1d / EBavgOcc1d), 2);
0513       EBlossMap2d[i - 1] = (EBloss_);
0514     }
0515
0516     TH2F* hEEmRecoMap2d(sAEReco.getME(0)->getTH2F());
0517     TH2F* hEEpRecoMap2d(sAEReco.getME(2)->getTH2F());
0518
0519     for (int i = 1; i < nEETowersPad - 1; i++) {
0520       for (int j = 0; j < nEEPhiTowers; j++) {
0521         int k = (i * nEETowersPad) + j + 1;
0522         int bin_ = hEEmRecoMap2d->GetBin(j + 1, i);
0523
0524         EEmrecoOcc1d[j] = eemRecoOccMap1dPad[k];
0525         EEminputOcc1d[j] = eemOccMap1dCumulPad[k];
0526         EEmavgOcc1d[j] = EEmavgOcc[k];
0527         double EEmcontent = eemRecoOccMap1dPad[k];
0528         hEEmRecoMap2d->SetBinContent(bin_, EEmcontent);
0529
0530         EEprecoOcc1d[j] = eepRecoOccMap1dPad[k];
0531         EEpinputOcc1d[j] = eepOccMap1dCumulPad[k];
0532         EEpavgOcc1d[j] = EEpavgOcc[k];
0533         double EEpcontent = eepRecoOccMap1dPad[k];
0534         hEEpRecoMap2d->SetBinContent(bin_, EEpcontent);
0535       }
0536       //Calculate the MSE loss = (output-input)^2, with avg response correction
0537       EEmloss_ = std::pow((EEmrecoOcc1d / EEmavgOcc1d - EEminputOcc1d / EEmavgOcc1d), 2);
0538       EEmlossMap2d[i - 1] = (EEmloss_);
0539
0540       EEploss_ = std::pow((EEprecoOcc1d / EEpavgOcc1d - EEpinputOcc1d / EEpavgOcc1d), 2);
0541       EEplossMap2d[i - 1] = (EEploss_);
0542     }
0543
0544     //Store each loss map from the output in the queue
0545     EBlossMap2dQ.push_back(EBlossMap2d);
0546     EEmlossMap2dQ.push_back(EEmlossMap2d);
0547     EEplossMap2dQ.push_back(EEplossMap2d);
0548
0549     //Keep exactly nLSloss loss maps to multiply
0550     if (EBlossMap2dQ.size() > nLSloss) {
0551       EBlossMap2dQ.pop_front();
0552       EEmlossMap2dQ.pop_front();
0553       EEplossMap2dQ.pop_front();
0554     }
0555     if (EBlossMap2dQ.size() < nLSloss) {  //Exit if there are not nLSloss loss maps
0556       return;
0557     }
0558
0559     //To hold the final multiplied loss
0560     std::valarray<std::valarray<float>> EBlossMap2dMult(std::valarray<float>(1., nEBPhiTowers), nEBEtaTowers);
0561     std::valarray<std::valarray<float>> EEmlossMap2dMult(std::valarray<float>(1., nEEPhiTowers), nEEEtaTowers);
0562     std::valarray<std::valarray<float>> EEplossMap2dMult(std::valarray<float>(1., nEEPhiTowers), nEEEtaTowers);
0563
0564     //Multiply together the last nLSloss loss maps
0565     //So that real anomalies which persist with time are enhanced and fluctuations are suppressed.
0566     for (size_t i = 0; i < EBlossMap2dQ.size(); i++) {
0567       EBlossMap2dMult *= EBlossMap2dQ[i];
0568       EEmlossMap2dMult *= EEmlossMap2dQ[i];
0569       EEplossMap2dMult *= EEplossMap2dQ[i];
0570     }
0571
0572     //Fill the AELoss ME with the values of this time multiplied loss map
0573     //MESet const& sAELoss(sources_.at("AELoss"));
0574     MESet& sAELoss(sources_.at("AELoss"));
0575
0576     TH2F* hEBLossMap2dMult(sAELoss.getME(1)->getTH2F());
0577
0578     for (int i = 0; i < hEBLossMap2dMult->GetNbinsY(); i++) {
0579       for (int j = 0; j < hEBLossMap2dMult->GetNbinsX(); j++) {
0580         int bin_ = hEBLossMap2dMult->GetBin(j + 1, i + 1);
0581         double content = EBlossMap2dMult[i][j];
0582         hEBLossMap2dMult->SetBinContent(bin_, content);
0583       }
0584     }
0585
0586     TH2F* hEEmLossMap2dMult(sAELoss.getME(0)->getTH2F());
0587     TH2F* hEEpLossMap2dMult(sAELoss.getME(2)->getTH2F());
0588
0589     for (int i = 0; i < hEEmLossMap2dMult->GetNbinsY(); i++) {
0590       for (int j = 0; j < hEEmLossMap2dMult->GetNbinsX(); j++) {
0591         int bin_ = hEEmLossMap2dMult->GetBin(j + 1, i + 1);
0592
0593         double EEmcontent = EEmlossMap2dMult[i][j];
0594         hEEmLossMap2dMult->SetBinContent(bin_, EEmcontent);
0595
0596         double EEpcontent = EEplossMap2dMult[i][j];
0597         hEEpLossMap2dMult->SetBinContent(bin_, EEpcontent);
0598       }
0599     }
0600
0601     ///////////////////// ML Quality Summary /////////////////////
0602     //Apply the quality threshold on the time multiplied loss map stored in the ME AELoss
0603     //If anomalous, the tower entry will have a large loss value. If good, the value will be close to zero.
0604
0605     MESet& meBadTowerCount(sources_.at("BadTowerCount"));
0606     MESet& meBadTowerCountNorm(sources_.at("BadTowerCountNorm"));
0607     MESet& meTrendMLBadTower(MEs_.at("TrendMLBadTower"));
0608
0609     LScount++;
0610
0611     MESet::const_iterator dAEnd(sAELoss.end(GetElectronicsMap()));
0612     for (MESet::const_iterator dItr(sAELoss.beginChannel(GetElectronicsMap())); dItr != dAEnd;
0613          dItr.toNextChannel(GetElectronicsMap())) {
0614       DetId id(dItr->getId());
0615
0616       bool doMaskML(meMLQualitySummary.maskMatches(id, mask, statusManager_, GetTrigTowerMap()));
0617
0618       float entries(dItr->getBinContent());
0619
0620       int quality(doMaskML ? kMGood : kGood);
0621       float MLThreshold;
0622
0623       if (id.subdetId() == EcalEndcap) {
0624         EEDetId eeid(id);
0625         if (eeid.zside() > 0)
0626           MLThreshold = EEpThreshold_;
0627         else
0628           MLThreshold = EEmThreshold_;
0629       } else {
0630         MLThreshold = EBThreshold_;
0631       }
0632
0633       //If a trigger tower entry is greater than the ML threshold, set it to Bad quality, otherwise Good.
0634       if (entries > MLThreshold) {
0635         quality = doMaskML ? kMBad : kBad;
0636         meBadTowerCount.fill(getEcalDQMSetupObjects(), id);
0637         if (id.subdetId() == EcalEndcap)
0638           nbadtowerEE++;
0639         else
0640           nbadtowerEB++;
0641       }
0642       //Fill the quality summary with the quality of the given tower id.
0643       meMLQualitySummary.setBinContent(getEcalDQMSetupObjects(), id, double(quality));
0644
0645       double badtowcount(meBadTowerCount.getBinContent(getEcalDQMSetupObjects(), id));
0646       meBadTowerCountNorm.setBinContent(getEcalDQMSetupObjects(), id, double(badtowcount / LScount));
0647     }  // ML Quality Summary
0648
0649     meTrendMLBadTower.fill(getEcalDQMSetupObjects(), EcalBarrel, double(timestamp_.iLumi), double(nbadtowerEB));
0650     meTrendMLBadTower.fill(getEcalDQMSetupObjects(), EcalEndcap, double(timestamp_.iLumi), double(nbadtowerEE));
0651
0652   }  // producePlots()
0653
0654   DEFINE_ECALDQM_WORKER(MLClient);
0655 }  // namespace ecaldqm