SiPixel/plugins/SiPixelGainCalibScaler.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    Tools/SiPixelGainCalibScaler
0004 // Class:      SiPixelGainCalibScaler
0005 //
0006 /**\class SiPixelGainCalibScaler SiPixelGainCalibScaler.cc Tools/SiPixelGainCalibScaler/plugins/SiPixelGainCalibScaler.cc
0007
0008  Description: Scales Pixel Gain Payloads by applying the VCal offset and slopes.
0009
0010  Implementation:
0011      Makes use of trick to loop over all IOVs in a tag by running on all the runs with EmptySource and just access DB once the IOV has changed via ESWatcher mechanism
0012 */
0013 //
0014 // Original Author:  Marco Musich
0015 //         Created:  Thu, 16 Jul 2020 10:36:21 GMT
0016 //
0017 //
0018
0019 // system include files
0020 #include <memory>
0021
0022 // user include files
0023 #include "CalibTracker/StandaloneTrackerTopology/interface/StandaloneTrackerTopology.h"
0024 #include "CondCore/DBOutputService/interface/PoolDBOutputService.h"
0025 #include "CondFormats/DataRecord/interface/SiPixelGainCalibrationForHLTRcd.h"
0026 #include "CondFormats/DataRecord/interface/SiPixelGainCalibrationOfflineRcd.h"
0027 #include "CondFormats/SiPixelObjects/interface/SiPixelGainCalibrationOffline.h"
0028 #include "CondFormats/SiPixelObjects/interface/SiPixelGainCalibrationForHLT.h"
0029 #include "DataFormats/DetId/interface/DetId.h"
0030 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0031 #include "DataFormats/TrackReco/interface/Track.h"
0032 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0033 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0034 #include "FWCore/Framework/interface/ESWatcher.h"
0035 #include "FWCore/Framework/interface/Event.h"
0036 #include "FWCore/Framework/interface/Frameworkfwd.h"
0037 #include "FWCore/Framework/interface/MakerMacros.h"
0038 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0039 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0040 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0041 #include "FWCore/ServiceRegistry/interface/Service.h"
0042 #include "FWCore/Utilities/interface/InputTag.h"
0043 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0044 #include "Geometry/CommonTopologies/interface/PixelTopology.h"
0045 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0046 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0047 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0048 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0049
0050 //
0051 // class declaration
0052 //
0053
0054 // If the analyzer does not use TFileService, please remove
0055 // the template argument to the base class so the class inherits
0056 // from  edm::one::EDAnalyzer<>
0057 // This will improve performance in multithreaded jobs.
0058
0059 namespace gainScale {
0060   struct VCalInfo {
0061   private:
0062     double m_conversionFactor;
0063     double m_conversionFactorL1;
0064     double m_offset;
0065     double m_offsetL1;
0066
0067   public:
0068     // default constructor
0069     VCalInfo() : m_conversionFactor(0.), m_conversionFactorL1(0.), m_offset(0.), m_offsetL1(0.) {}
0070
0071     // initialize
0072     void init(double conversionFactor, double conversionFactorL1, double offset, double offsetL1) {
0073       m_conversionFactor = conversionFactor;
0074       m_conversionFactorL1 = conversionFactorL1;
0075       m_offset = offset;
0076       m_offsetL1 = offsetL1;
0077     }
0078
0079     void printAllInfo() {
0080       edm::LogVerbatim("SiPixelGainCalibScaler") << " conversion factor      : " << m_conversionFactor << "\n"
0081                                                  << " conversion factor (L1) : " << m_conversionFactorL1 << "\n"
0082                                                  << " offset                 : " << m_offset << "\n"
0083                                                  << " offset            (L1) : " << m_offsetL1 << "\n";
0084     }
0085
0086     double getConversionFactor() { return m_conversionFactor; }
0087     double getConversionFactorL1() { return m_conversionFactorL1; }
0088     double getOffset() { return m_offset; }
0089     double getOffsetL1() { return m_offsetL1; }
0090     virtual ~VCalInfo() {}
0091   };
0092 }  // namespace gainScale
0093
0094 class SiPixelGainCalibScaler : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0095 public:
0096   explicit SiPixelGainCalibScaler(const edm::ParameterSet&);
0097   ~SiPixelGainCalibScaler() override;
0098
0099   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0100
0101 private:
0102   void beginJob() override;
0103   void analyze(const edm::Event&, const edm::EventSetup&) override;
0104   void endJob() override;
0105   template <class tokenType, class PayloadType>
0106   void computeAndStorePalyoads(const edm::EventSetup& iSetup, const tokenType& token);
0107
0108   // ----------member data ---------------------------
0109   const std::string recordName_;
0110   const bool isForHLT_;
0111   const bool verbose_;
0112   const std::vector<edm::ParameterSet> m_parameters;
0113
0114   gainScale::VCalInfo phase0VCal;
0115   gainScale::VCalInfo phase1VCal;
0116
0117   edm::ESGetToken<SiPixelGainCalibrationForHLT, SiPixelGainCalibrationForHLTRcd> gainHLTCalibToken_;
0118   edm::ESGetToken<SiPixelGainCalibrationOffline, SiPixelGainCalibrationOfflineRcd> gainOfflineCalibToken_;
0119   edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> tkGeomToken_;
0120   edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> tkTopoToken_;
0121
0122   edm::ESWatcher<SiPixelGainCalibrationForHLTRcd> pixelHLTGainWatcher_;
0123   edm::ESWatcher<SiPixelGainCalibrationOfflineRcd> pixelOfflineGainWatcher_;
0124 };
0125
0126 //
0127 // constructors and destructor
0128 //
0129 SiPixelGainCalibScaler::SiPixelGainCalibScaler(const edm::ParameterSet& iConfig)
0130     : recordName_(iConfig.getParameter<std::string>("record")),
0131       isForHLT_(iConfig.getParameter<bool>("isForHLT")),
0132       verbose_(iConfig.getUntrackedParameter<bool>("verbose", false)),
0133       m_parameters(iConfig.getParameter<std::vector<edm::ParameterSet> >("parameters")) {
0134   gainHLTCalibToken_ = esConsumes();
0135   gainOfflineCalibToken_ = esConsumes();
0136   tkGeomToken_ = esConsumes();
0137   tkTopoToken_ = esConsumes();
0138
0139   for (auto& thePSet : m_parameters) {
0140     const unsigned int phase(thePSet.getParameter<unsigned int>("phase"));
0141     switch (phase) {
0142       case 0: {
0143         phase0VCal.init(thePSet.getParameter<double>("conversionFactor"),
0144                         thePSet.getParameter<double>("conversionFactorL1"),
0145                         thePSet.getParameter<double>("offset"),
0146                         thePSet.getParameter<double>("offsetL1"));
0147         break;
0148       }
0149       case 1: {
0150         phase1VCal.init(thePSet.getParameter<double>("conversionFactor"),
0151                         thePSet.getParameter<double>("conversionFactorL1"),
0152                         thePSet.getParameter<double>("offset"),
0153                         thePSet.getParameter<double>("offsetL1"));
0154         break;
0155       }
0156       default:
0157         throw cms::Exception("LogicError") << "Unrecongnized phase: " << phase << ". Exiting!";
0158     }
0159   }
0160 }
0161
0162 SiPixelGainCalibScaler::~SiPixelGainCalibScaler() {}
0163
0164 //
0165 // member functions
0166 //
0167
0168 // ------------ method called for each event  ------------
0169 void SiPixelGainCalibScaler::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0170   using namespace edm;
0171
0172   int run = iEvent.id().run();
0173   bool hasPixelHLTGainIOV = pixelHLTGainWatcher_.check(iSetup);
0174   bool hasPixelOfflineGainIOV = pixelOfflineGainWatcher_.check(iSetup);
0175
0176   if ((hasPixelHLTGainIOV && isForHLT_) || (hasPixelOfflineGainIOV && !isForHLT_)) {
0177     edm::LogPrint("SiPixelGainCalibScaler") << " Pixel Gains have a new IOV for run: " << run << std::endl;
0178
0179     if (isForHLT_) {
0180       computeAndStorePalyoads<edm::ESGetToken<SiPixelGainCalibrationForHLT, SiPixelGainCalibrationForHLTRcd>,
0181                               SiPixelGainCalibrationForHLT>(iSetup, gainHLTCalibToken_);
0182     } else {
0183       computeAndStorePalyoads<edm::ESGetToken<SiPixelGainCalibrationOffline, SiPixelGainCalibrationOfflineRcd>,
0184                               SiPixelGainCalibrationOffline>(iSetup, gainOfflineCalibToken_);
0185     }
0186   }  // if new IOV
0187 }
0188
0189 // ------------ template method to construct the payloads  ------------
0190 template <class tokenType, class PayloadType>
0191 void SiPixelGainCalibScaler::computeAndStorePalyoads(const edm::EventSetup& iSetup, const tokenType& token) {
0192   gainScale::VCalInfo myVCalInfo;
0193
0194   //=======================================================
0195   // Retrieve geometry information
0196   //=======================================================
0197   const TrackerGeometry* pDD = &iSetup.getData(tkGeomToken_);
0198   edm::LogInfo("SiPixelGainCalibScaler") << "There are: " << pDD->dets().size() << " detectors";
0199
0200   // switch on the phase1
0201   if ((pDD->isThere(GeomDetEnumerators::P1PXB)) || (pDD->isThere(GeomDetEnumerators::P1PXEC))) {
0202     myVCalInfo = phase1VCal;
0203     edm::LogInfo("SiPixelGainCalibScaler") << " ==> This is a phase1 IOV";
0204   } else {
0205     myVCalInfo = phase0VCal;
0206     edm::LogInfo("SiPixelGainCalibScaler") << " ==> This is a phase0 IOV";
0207   }
0208
0209   myVCalInfo.printAllInfo();
0210
0211   // if need the ESHandle to check if the SetupData was there or not
0212   auto payload = &iSetup.getData(token);
0213   std::vector<uint32_t> detids;
0214   payload->getDetIds(detids);
0215
0216   float mingain = payload->getGainLow();
0217   float maxgain = (payload->getGainHigh()) * myVCalInfo.getConversionFactorL1();
0218   float minped = payload->getPedLow();
0219   float maxped = payload->getPedHigh() * 1.10;
0220
0221   PayloadType SiPixelGainCalibration_(minped, maxped, mingain, maxgain);
0222
0223   //Retrieve tracker topology from geometry
0224   const TrackerTopology* tTopo = &iSetup.getData(tkTopoToken_);
0225
0226   // possible to load it not from EventSetup
0227   //const char* path_toTopologyXML = "Geometry/TrackerCommonData/data/PhaseI/trackerParameters.xml";
0228   //TrackerTopology tTopo = StandaloneTrackerTopology::fromTrackerParametersXMLFile(edm::FileInPath(path_toTopologyXML).fullPath());
0229
0230   for (const auto& d : detids) {
0231     bool isLayer1 = false;
0232     int subid = DetId(d).subdetId();
0233     if (subid == PixelSubdetector::PixelBarrel) {
0234       auto layer = tTopo->pxbLayer(DetId(d));
0235       if (layer == 1) {
0236         isLayer1 = true;
0237       }
0238     }
0239
0240     std::vector<char> theSiPixelGainCalibration;
0241
0242     auto range = payload->getRange(d);
0243     int numberOfRowsToAverageOver = payload->getNumberOfRowsToAverageOver();
0244     int ncols = payload->getNCols(d);
0245     int nRocsInRow = (range.second - range.first) / ncols / numberOfRowsToAverageOver;
0246     unsigned int nRowsForHLT = 1;
0247     int nrows = std::max((payload->getNumberOfRowsToAverageOver() * nRocsInRow),
0248                          nRowsForHLT);  // dirty trick to make it work for the HLT payload
0249
0250     auto rangeAndCol = payload->getRangeAndNCols(d);
0251     bool isDeadColumn;
0252     bool isNoisyColumn;
0253
0254     if (verbose_) {
0255       edm::LogVerbatim("SiPixelGainCalibScaler")
0256           << "NCOLS: " << payload->getNCols(d) << " " << rangeAndCol.second << " NROWS:" << nrows
0257           << ", RANGES: " << rangeAndCol.first.second - rangeAndCol.first.first
0258           << ", Ratio: " << float(rangeAndCol.first.second - rangeAndCol.first.first) / rangeAndCol.second << std::endl;
0259     }
0260
0261     for (int col = 0; col < ncols; col++) {
0262       for (int row = 0; row < nrows; row++) {
0263         float gain = payload->getGain(col, row, rangeAndCol.first, rangeAndCol.second, isDeadColumn, isNoisyColumn);
0264         float ped = payload->getPed(col, row, rangeAndCol.first, rangeAndCol.second, isDeadColumn, isNoisyColumn);
0265
0266         if (verbose_)
0267           edm::LogInfo("SiPixelGainCalibScaler") << "pre-change gain: " << gain << " pede:" << ped << std::endl;
0268
0269         //
0270         // From here https://github.com/cms-sw/cmssw/blob/master/CalibTracker/SiPixelESProducers/src/SiPixelGainCalibrationForHLTService.cc#L20-L47
0271         //
0272         // vcal = ADC * DBgain - DBped * DBgain
0273         // electrons = vcal * conversionFactor + offset
0274         //
0275         // follows:
0276         // electrons = (ADC*DBgain � DBped*DBgain)*conversionFactor + offset
0277         // electrons = ADC*conversionFactor*DBgain - conversionFactor*DBped*DBgain + offset
0278         //
0279         // this should equal the new equation:
0280         //
0281         // electrons = ADC*DBgain' - DBPed' * DBgain'
0282         //
0283         // So equating piece by piece:
0284         //
0285         // DBgain' = conversionFactor*DBgain
0286         // DBped' = (conversionFactor*DBped*Dbgain � offset)/(conversionFactor*DBgain)
0287         //        = DBped - offset/DBgain'
0288         //
0289
0290         if (isLayer1) {
0291           gain = gain * myVCalInfo.getConversionFactorL1();
0292           ped = ped - myVCalInfo.getOffsetL1() / gain;
0293         } else {
0294           gain = gain * myVCalInfo.getConversionFactor();
0295           ped = ped - myVCalInfo.getOffset() / gain;
0296         }
0297
0298         if (verbose_)
0299           edm::LogInfo("SiPixelGainCalibScaler") << "post-change gain: " << gain << " pede:" << ped << std::endl;
0300
0301         if constexpr (std::is_same_v<PayloadType, SiPixelGainCalibrationForHLT>) {
0302           SiPixelGainCalibration_.setData(ped, gain, theSiPixelGainCalibration, false, false);
0303         } else {
0304           SiPixelGainCalibration_.setDataPedestal(ped, theSiPixelGainCalibration);
0305           if ((row + 1) % numberOfRowsToAverageOver == 0) {  // fill the column average after every ROC!
0306             SiPixelGainCalibration_.setDataGain(gain, numberOfRowsToAverageOver, theSiPixelGainCalibration);
0307           }
0308         }
0309       }  // loop on rows
0310     }  // loop on columns
0311
0312     typename PayloadType::Range outrange(theSiPixelGainCalibration.begin(), theSiPixelGainCalibration.end());
0313     if (!SiPixelGainCalibration_.put(d, outrange, ncols))
0314       edm::LogError("SiPixelGainCalibScaler") << "[SiPixelGainCalibScaler::analyze] detid already exists" << std::endl;
0315   }  // loop on DetIds
0316
0317   // Write into DB
0318   edm::LogInfo(" --- writing to DB!");
0319   edm::Service<cond::service::PoolDBOutputService> mydbservice;
0320   if (!mydbservice.isAvailable()) {
0321     edm::LogError("db service unavailable");
0322     return;
0323   } else {
0324     edm::LogInfo("DB service OK");
0325   }
0326
0327   try {
0328     if (mydbservice->isNewTagRequest(recordName_)) {
0329       mydbservice->createOneIOV<PayloadType>(SiPixelGainCalibration_, mydbservice->beginOfTime(), recordName_);
0330     } else {
0331       mydbservice->appendOneIOV<PayloadType>(SiPixelGainCalibration_, mydbservice->currentTime(), recordName_);
0332     }
0333     edm::LogInfo(" --- all OK");
0334   } catch (const cond::Exception& er) {
0335     edm::LogError("SiPixelGainCalibScaler") << er.what() << std::endl;
0336   } catch (const std::exception& er) {
0337     edm::LogError("SiPixelGainCalibScaler") << "caught std::exception " << er.what() << std::endl;
0338   } catch (...) {
0339     edm::LogError("SiPixelGainCalibScaler") << "Funny error" << std::endl;
0340   }
0341 }
0342
0343 // ------------ method called once each job just before starting event loop  ------------
0344 void SiPixelGainCalibScaler::beginJob() {}
0345
0346 // ------------ method called once each job just after ending the event loop  ------------
0347 void SiPixelGainCalibScaler::endJob() {}
0348
0349 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0350 void SiPixelGainCalibScaler::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0351   edm::ParameterSetDescription desc;
0352   desc.add<std::string>("record", "SiPixelGainCalibrationForHLTRcd");
0353   desc.add<bool>("isForHLT", true);
0354
0355   edm::ParameterSetDescription vcalInfos;
0356   vcalInfos.add<unsigned int>("phase");
0357   vcalInfos.add<double>("conversionFactor");
0358   vcalInfos.add<double>("conversionFactorL1");
0359   vcalInfos.add<double>("offset");
0360   vcalInfos.add<double>("offsetL1");
0361
0362   std::vector<edm::ParameterSet> tmp;
0363   tmp.reserve(2);
0364   {
0365     edm::ParameterSet phase0VCal;
0366     phase0VCal.addParameter<unsigned int>("phase", 0);
0367     phase0VCal.addParameter<double>("conversionFactor", 65.);
0368     phase0VCal.addParameter<double>("conversionFactorL1", 65.);
0369     phase0VCal.addParameter<double>("offset", -414.);
0370     phase0VCal.addParameter<double>("offsetL1", -414.);
0371     tmp.push_back(phase0VCal);
0372   }
0373   {
0374     edm::ParameterSet phase1VCal;
0375     phase1VCal.addParameter<unsigned int>("phase", 1);
0376     phase1VCal.addParameter<double>("conversionFactor", 47.);
0377     phase1VCal.addParameter<double>("conversionFactorL1", 50.);
0378     phase1VCal.addParameter<double>("offset", -60.);
0379     phase1VCal.addParameter<double>("offsetL1", -670.);
0380     tmp.push_back(phase1VCal);
0381   }
0382   desc.addVPSet("parameters", vcalInfos, tmp);
0383
0384   desc.addUntracked<bool>("verbose", false);
0385
0386   descriptions.add("siPixelGainCalibScaler", desc);
0387 }
0388
0389 //define this as a plug-in
0390 DEFINE_FWK_MODULE(SiPixelGainCalibScaler);