plugins/alpaka/HGCalRecHitCalibrationESProducer.cc

0001 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0002 #include "FWCore/ParameterSet/interface/FileInPath.h"
0003 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0004
0005 #include "HeterogeneousCore/AlpakaCore/interface/alpaka/ESGetToken.h"
0006 #include "HeterogeneousCore/AlpakaCore/interface/alpaka/ESProducer.h"
0007 #include "HeterogeneousCore/AlpakaCore/interface/alpaka/ModuleFactory.h"
0008 #include "HeterogeneousCore/AlpakaInterface/interface/config.h"
0009 #include "HeterogeneousCore/AlpakaInterface/interface/host.h"
0010 #include "HeterogeneousCore/AlpakaInterface/interface/memory.h"
0011
0012 #include "DataFormats/HGCalDigi/interface/HGCalElectronicsId.h"
0013 #include "CondFormats/HGCalObjects/interface/HGCalMappingModuleIndexer.h"
0014 #include "CondFormats/HGCalObjects/interface/HGCalCalibrationParameterHost.h"
0015 #include "CondFormats/HGCalObjects/interface/alpaka/HGCalCalibrationParameterDevice.h"
0016 #include "CondFormats/DataRecord/interface/HGCalElectronicsMappingRcd.h"
0017 #include "CondFormats/DataRecord/interface/HGCalModuleConfigurationRcd.h"  // depends on HGCalElectronicsMappingRcd
0018 #include "RecoLocalCalo/HGCalRecAlgos/interface/HGCalESProducerTools.h"    // for json, search_modkey
0019
0020 #include <string>
0021 //#include <iostream>
0022 //#include <sstream>
0023 #include <fstream>  // needed to read json file with std::ifstream
0024
0025 namespace ALPAKA_ACCELERATOR_NAMESPACE {
0026
0027   namespace hgcalrechit {
0028
0029     class HGCalCalibrationESProducer : public ESProducer {
0030     public:
0031       HGCalCalibrationESProducer(const edm::ParameterSet& iConfig)
0032           : ESProducer(iConfig), filename_(iConfig.getParameter<edm::FileInPath>("filename")) {
0033         auto cc = setWhatProduced(this);
0034         indexToken_ = cc.consumes(iConfig.getParameter<edm::ESInputTag>("indexSource"));
0035       }
0036
0037       static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0038         edm::ParameterSetDescription desc;
0039         desc.add<edm::FileInPath>("filename")->setComment("Path to JSON file with calibration parameters");
0040         desc.add<edm::ESInputTag>("indexSource", edm::ESInputTag(""))
0041             ->setComment("Label for module indexer to set SoA size");
0042         descriptions.addWithDefaultLabel(desc);
0043       }
0044
0045       // @short fill SoA column with data from vector for any type with some offset
0046       template <typename T>
0047       static void fill_SoA_column(
0048           T* column_SoA, const std::vector<T>& values, const int offset, const int nrows, int arr_offset = 0) {
0049         const int nrows_vals = values.size();
0050         if (arr_offset < 0) {
0051           arr_offset = 0;
0052           if (nrows_vals != arr_offset + nrows) {
0053             throw edm::Exception(edm::errors::LogicError, "HGCalCalibrationESProducer")
0054                 << " Expected " << nrows << " rows, but got " << nrows_vals << "!";
0055           }
0056         } else if (nrows_vals < arr_offset + nrows) {
0057           throw edm::Exception(edm::errors::LogicError, "HGCalCalibrationESProducer")
0058               << " Tried to copy " << nrows << " rows to SoA with offset " << arr_offset << ", but only have "
0059               << nrows_vals << " values in JSON!";
0060         }
0061         auto begin = values.begin() + arr_offset;
0062         auto end = (begin + nrows > values.end()) ? values.end() : begin + nrows;
0063         std::copy(begin, end, &column_SoA[offset]);
0064       }
0065
0066       // @short fill full SoA column with data from vector for any type
0067       template <typename T, typename P>
0068       static void fill_SoA_eigen_row(P& soa, const std::vector<std::vector<T>>& values, const size_t row) {
0069         if (row >= values.size())
0070           throw edm::Exception(edm::errors::LogicError, "HGCalCalibrationESProducer")
0071               << " Tried to copy row " << row << " to SoA, but only have " << values.size() << " values in JSON!";
0072         if (!values.empty() && int(values[row].size()) != soa.size())
0073           throw edm::Exception(edm::errors::LogicError, "HGCalCalibrationESProducer")
0074               << " Expected " << soa.size() << " elements in Eigen vector, but got " << values[row].size() << "!";
0075         for (int i = 0; i < soa.size(); i++)
0076           soa(i) = values[row][i];
0077       }
0078
0079       // @short create the ESProducer product: a SoA with channel-level calibration constants
0080       std::optional<hgcalrechit::HGCalCalibParamHost> produce(const HGCalModuleConfigurationRcd& iRecord) {
0081         auto const& moduleMap = iRecord.get(indexToken_);
0082         edm::LogInfo("HGCalCalibrationESProducer") << "produce: filename=" << filename_.fullPath().c_str();
0083
0084         // load dense indexing
0085         const uint32_t nchans = moduleMap.getMaxDataSize();  // channel-level size
0086         hgcalrechit::HGCalCalibParamHost product(nchans, cms::alpakatools::host());
0087
0088         // load calib parameters from JSON
0089         std::ifstream infile(filename_.fullPath().c_str());
0090         json calib_data = json::parse(infile, nullptr, true, /*ignore_comments*/ true);
0091         for (const auto& it : moduleMap.getTypecodeMap()) {  // loop over all module typecodes
0092           std::string const& module = it.first;              // module typecode, e.g. "ML-F3PT-TX-0003"
0093
0094           // retrieve matching key (glob patterns allowed)
0095           const auto modkey = hgcal::search_modkey(module, calib_data, filename_.fullPath());
0096           auto calib_data_ = calib_data[modkey];
0097
0098           // get dimensions
0099           const auto firstkey = calib_data_.begin().key();
0100           const uint32_t offset = moduleMap.getIndexForModuleData(module);  // first channel index
0101           const uint32_t nchans = moduleMap.getNumChannels(module);         // number of channels in mapper
0102           uint32_t nrows = calib_data_[firstkey].size();                    // number of channels in JSON
0103
0104           // check number of channels & ROCs make sense
0105           if (nrows % 37 != 0) {
0106             edm::LogWarning("HGCalCalibrationESProducer")
0107                 << "nchannels=" << nrows << ", which is not divisible by 37 (#channels per e-Rx)!";
0108           }
0109           if (nchans != nrows) {
0110             edm::LogWarning("HGCalCalibrationESProducer")
0111                 << "nchannels does not match between module indexer ('" << module << "'," << nchans << ") and JSON ('"
0112                 << modkey << "'," << nrows << "')!";
0113             nrows = std::min(nrows, nchans);  // take smallest to avoid overlap
0114           }
0115
0116           // fill calibration parameters for ADC, CM, TOT, MIPS scale, ...
0117           fill_SoA_column<float>(product.view().ADC_ped(), calib_data_["ADC_ped"], offset, nrows);
0118           fill_SoA_column<float>(product.view().Noise(), calib_data_["Noise"], offset, nrows);
0119           fill_SoA_column<float>(product.view().CM_slope(), calib_data_["CM_slope"], offset, nrows);
0120           fill_SoA_column<float>(product.view().CM_ped(), calib_data_["CM_ped"], offset, nrows);
0121           fill_SoA_column<float>(product.view().BXm1_slope(), calib_data_["BXm1_slope"], offset, nrows);
0122           fill_SoA_column<float>(product.view().TOTtoADC(), calib_data_["TOTtoADC"], offset, nrows);
0123           fill_SoA_column<float>(product.view().TOT_ped(), calib_data_["TOT_ped"], offset, nrows);
0124           fill_SoA_column<float>(product.view().TOT_lin(), calib_data_["TOT_lin"], offset, nrows);
0125           fill_SoA_column<float>(product.view().TOT_P0(), calib_data_["TOT_P0"], offset, nrows);
0126           fill_SoA_column<float>(product.view().TOT_P1(), calib_data_["TOT_P1"], offset, nrows);
0127           fill_SoA_column<float>(product.view().TOT_P2(), calib_data_["TOT_P2"], offset, nrows);
0128           fill_SoA_column<float>(product.view().MIPS_scale(), calib_data_["MIPS_scale"], offset, nrows);
0129           fill_SoA_column<unsigned char>(product.view().valid(), calib_data_["Valid"], offset, nrows);
0130
0131           // default parameters for fine-grain TDC, coarse-grain TDC and time-walk corrections that allow passthrough
0132           // (for backwards compatibility of older calibration JSONs)
0133           if (calib_data_.find("TOA_CTDC") == calib_data_.end())
0134             calib_data_["TOA_CTDC"] = std::vector<std::vector<float>>(nrows, std::vector<float>(32, 0.));
0135           if (calib_data_.find("TOA_FTDC") == calib_data_.end())
0136             calib_data_["TOA_FTDC"] = std::vector<std::vector<float>>(nrows, std::vector<float>(8, 0.));
0137           if (calib_data_.find("TOA_TW") == calib_data_.end())
0138             calib_data_["TOA_TW"] = std::vector<std::vector<float>>(nrows, std::vector<float>(3, 0.));
0139
0140           // fill vectors for ToA correction parameters
0141           for (size_t n = 0; n < nrows; n++) {
0142             auto vi = product.view()[offset + n];
0143             fill_SoA_eigen_row<float>(vi.TOA_CTDC(), calib_data_["TOA_CTDC"], n);
0144             fill_SoA_eigen_row<float>(vi.TOA_FTDC(), calib_data_["TOA_FTDC"], n);
0145             fill_SoA_eigen_row<float>(vi.TOA_TW(), calib_data_["TOA_TW"], n);
0146           }
0147         }
0148
0149         return product;
0150       }  // end of produce()
0151
0152     private:
0153       edm::ESGetToken<HGCalMappingModuleIndexer, HGCalElectronicsMappingRcd> indexToken_;
0154       const edm::FileInPath filename_;
0155     };
0156
0157   }  // namespace hgcalrechit
0158
0159 }  // namespace ALPAKA_ACCELERATOR_NAMESPACE
0160
0161 DEFINE_FWK_EVENTSETUP_ALPAKA_MODULE(hgcalrechit::HGCalCalibrationESProducer);