diff --git a/CondTools/SiPixel/test/SiPixelGainCalibScaler.cc b/CondTools/SiPixel/test/SiPixelGainCalibScaler.cc new file mode 100644 index 0000000000000..6b2388d85887c --- /dev/null +++ b/CondTools/SiPixel/test/SiPixelGainCalibScaler.cc @@ -0,0 +1,389 @@ +// -*- C++ -*- +// +// Package: Tools/SiPixelGainCalibScaler +// Class: SiPixelGainCalibScaler +// +/**\class SiPixelGainCalibScaler SiPixelGainCalibScaler.cc Tools/SiPixelGainCalibScaler/plugins/SiPixelGainCalibScaler.cc + + Description: Scales Pixel Gain Payloads by applying the VCal offset and slopes. + + Implementation: + 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 +*/ +// +// Original Author: Marco Musich +// Created: Thu, 16 Jul 2020 10:36:21 GMT +// +// + +// system include files +#include + +// user include files +#include "CalibTracker/StandaloneTrackerTopology/interface/StandaloneTrackerTopology.h" +#include "CondCore/DBOutputService/interface/PoolDBOutputService.h" +#include "CondFormats/DataRecord/interface/SiPixelGainCalibrationForHLTRcd.h" +#include "CondFormats/DataRecord/interface/SiPixelGainCalibrationOfflineRcd.h" +#include "CondFormats/SiPixelObjects/interface/SiPixelGainCalibrationOffline.h" +#include "CondFormats/SiPixelObjects/interface/SiPixelGainCalibrationForHLT.h" +#include "DataFormats/DetId/interface/DetId.h" +#include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h" +#include "DataFormats/TrackReco/interface/Track.h" +#include "DataFormats/TrackReco/interface/TrackFwd.h" +#include "DataFormats/TrackerCommon/interface/TrackerTopology.h" +#include "FWCore/Framework/interface/ESWatcher.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/Framework/interface/one/EDAnalyzer.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/ServiceRegistry/interface/Service.h" +#include "FWCore/Utilities/interface/InputTag.h" +#include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h" +#include "Geometry/CommonTopologies/interface/PixelTopology.h" +#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h" +#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h" +#include "Geometry/Records/interface/TrackerTopologyRcd.h" +#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h" + +// +// class declaration +// + +// If the analyzer does not use TFileService, please remove +// the template argument to the base class so the class inherits +// from edm::one::EDAnalyzer<> +// This will improve performance in multithreaded jobs. + +namespace gainScale { + struct VCalInfo { + private: + double m_conversionFactor; + double m_conversionFactorL1; + double m_offset; + double m_offsetL1; + + public: + // default constructor + VCalInfo() : m_conversionFactor(0.), m_conversionFactorL1(0.), m_offset(0.), m_offsetL1(0.) {} + + // initialize + void init(double conversionFactor, double conversionFactorL1, double offset, double offsetL1) { + m_conversionFactor = conversionFactor; + m_conversionFactorL1 = conversionFactorL1; + m_offset = offset; + m_offsetL1 = offsetL1; + } + + void printAllInfo() { + edm::LogVerbatim("SiPixelGainCalibScaler") << " conversion factor : " << m_conversionFactor << "\n" + << " conversion factor (L1) : " << m_conversionFactorL1 << "\n" + << " offset : " << m_offset << "\n" + << " offset (L1) : " << m_offsetL1 << "\n"; + } + + double getConversionFactor() { return m_conversionFactor; } + double getConversionFactorL1() { return m_conversionFactorL1; } + double getOffset() { return m_offset; } + double getOffsetL1() { return m_offsetL1; } + virtual ~VCalInfo() {} + }; +} // namespace gainScale + +class SiPixelGainCalibScaler : public edm::one::EDAnalyzer { +public: + explicit SiPixelGainCalibScaler(const edm::ParameterSet&); + ~SiPixelGainCalibScaler() override; + + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); + +private: + void beginJob() override; + void analyze(const edm::Event&, const edm::EventSetup&) override; + void endJob() override; + template + void computeAndStorePalyoads(const edm::EventSetup& iSetup, const tokenType& token); + + // ----------member data --------------------------- + const std::string recordName_; + const bool isForHLT_; + const bool verbose_; + const std::vector m_parameters; + + gainScale::VCalInfo phase0VCal; + gainScale::VCalInfo phase1VCal; + + edm::ESGetToken gainHLTCalibToken_; + edm::ESGetToken gainOfflineCalibToken_; + + edm::ESWatcher pixelHLTGainWatcher_; + edm::ESWatcher pixelOfflineGainWatcher_; +}; + +// +// constructors and destructor +// +SiPixelGainCalibScaler::SiPixelGainCalibScaler(const edm::ParameterSet& iConfig) + : recordName_(iConfig.getParameter("record")), + isForHLT_(iConfig.getParameter("isForHLT")), + verbose_(iConfig.getUntrackedParameter("verbose", false)), + m_parameters(iConfig.getParameter >("parameters")) { + gainHLTCalibToken_ = esConsumes(); + gainOfflineCalibToken_ = esConsumes(); + + for (auto& thePSet : m_parameters) { + const unsigned int phase(thePSet.getParameter("phase")); + switch (phase) { + case 0: { + phase0VCal.init(thePSet.getParameter("conversionFactor"), + thePSet.getParameter("conversionFactorL1"), + thePSet.getParameter("offset"), + thePSet.getParameter("offsetL1")); + break; + } + case 1: { + phase1VCal.init(thePSet.getParameter("conversionFactor"), + thePSet.getParameter("conversionFactorL1"), + thePSet.getParameter("offset"), + thePSet.getParameter("offsetL1")); + break; + } + default: + throw cms::Exception("LogicError") << "Unrecongnized phase: " << phase << ". Exiting!"; + } + } +} + +SiPixelGainCalibScaler::~SiPixelGainCalibScaler() {} + +// +// member functions +// + +// ------------ method called for each event ------------ +void SiPixelGainCalibScaler::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) { + using namespace edm; + + int run = iEvent.id().run(); + bool hasPixelHLTGainIOV = pixelHLTGainWatcher_.check(iSetup); + bool hasPixelOfflineGainIOV = pixelOfflineGainWatcher_.check(iSetup); + + if ((hasPixelHLTGainIOV && isForHLT_) || (hasPixelOfflineGainIOV && !isForHLT_)) { + edm::LogPrint("SiPixelGainCalibScaler") << " Pixel Gains have a new IOV for run: " << run << std::endl; + + if (isForHLT_) { + computeAndStorePalyoads, + SiPixelGainCalibrationForHLT>(iSetup, gainHLTCalibToken_); + } else { + computeAndStorePalyoads, + SiPixelGainCalibrationOffline>(iSetup, gainOfflineCalibToken_); + } + } // if new IOV +} + +// ------------ template method to construct the payloads ------------ +template +void SiPixelGainCalibScaler::computeAndStorePalyoads(const edm::EventSetup& iSetup, const tokenType& token) { + gainScale::VCalInfo myVCalInfo; + + //======================================================= + // Retrieve geometry information + //======================================================= + edm::ESHandle pDD; + iSetup.get().get(pDD); + edm::LogInfo("SiPixelGainCalibScaler") << "There are: " << pDD->dets().size() << " detectors"; + + // switch on the phase1 + if ((pDD->isThere(GeomDetEnumerators::P1PXB)) || (pDD->isThere(GeomDetEnumerators::P1PXEC))) { + myVCalInfo = phase1VCal; + edm::LogInfo("SiPixelGainCalibScaler") << " ==> This is a phase1 IOV"; + } else { + myVCalInfo = phase0VCal; + edm::LogInfo("SiPixelGainCalibScaler") << " ==> This is a phase0 IOV"; + } + + myVCalInfo.printAllInfo(); + + // if need the ESHandle to check if the SetupData was there or not + auto payload = iSetup.getHandle(token); + std::vector detids; + payload->getDetIds(detids); + + float mingain = payload->getGainLow(); + float maxgain = (payload->getGainHigh()) * myVCalInfo.getConversionFactorL1(); + float minped = payload->getPedLow(); + float maxped = payload->getPedHigh() * 1.10; + + auto SiPixelGainCalibration_ = new PayloadType(minped, maxped, mingain, maxgain); + + //Retrieve tracker topology from geometry + edm::ESHandle tTopoHandle; + iSetup.get().get(tTopoHandle); + const TrackerTopology* tTopo = tTopoHandle.product(); + + //const char* path_toTopologyXML = "Geometry/TrackerCommonData/data/PhaseI/trackerParameters.xml"; + //TrackerTopology tTopo = StandaloneTrackerTopology::fromTrackerParametersXMLFile(edm::FileInPath(path_toTopologyXML).fullPath()); + + for (const auto& d : detids) { + bool isLayer1 = false; + int subid = DetId(d).subdetId(); + if (subid == PixelSubdetector::PixelBarrel) { + auto layer = tTopo->pxbLayer(DetId(d)); + if (layer == 1) { + isLayer1 = true; + } + } + + std::vector theSiPixelGainCalibration; + + auto range = payload->getRange(d); + int numberOfRowsToAverageOver = payload->getNumberOfRowsToAverageOver(); + int ncols = payload->getNCols(d); + int nRocsInRow = (range.second - range.first) / ncols / numberOfRowsToAverageOver; + unsigned int nRowsForHLT = 1; + int nrows = std::max((payload->getNumberOfRowsToAverageOver() * nRocsInRow), + nRowsForHLT); // dirty trick to make it work for the HLT payload + + auto rangeAndCol = payload->getRangeAndNCols(d); + bool isDeadColumn; + bool isNoisyColumn; + + if (verbose_) { + edm::LogVerbatim("SiPixelGainCalibScaler") + << "NCOLS: " << payload->getNCols(d) << " " << rangeAndCol.second << " NROWS:" << nrows + << ", RANGES: " << rangeAndCol.first.second - rangeAndCol.first.first + << ", Ratio: " << float(rangeAndCol.first.second - rangeAndCol.first.first) / rangeAndCol.second << std::endl; + } + + for (int col = 0; col < ncols; col++) { + for (int row = 0; row < nrows; row++) { + float gain = payload->getGain(col, row, rangeAndCol.first, rangeAndCol.second, isDeadColumn, isNoisyColumn); + float ped = payload->getPed(col, row, rangeAndCol.first, rangeAndCol.second, isDeadColumn, isNoisyColumn); + + if (verbose_) + edm::LogInfo("SiPixelGainCalibScaler") << "pre-change gain: " << gain << " pede:" << ped << std::endl; + + // + // From here https://github.com/cms-sw/cmssw/blob/master/CalibTracker/SiPixelESProducers/src/SiPixelGainCalibrationForHLTService.cc#L20-L47 + // + // vcal = ADC * DBgain - DBped * DBgain + // electrons = vcal * conversionFactor + offset + // + // follows: + // electrons = (ADC*DBgain � DBped*DBgain)*conversionFactor + offset + // electrons = ADC*conversionFactor*DBgain - conversionFactor*DBped*DBgain + offset + // + // this should equal the new equation: + // + // electrons = ADC*DBgain' - DBPed' * DBgain' + // + // So equating piece by piece: + // + // DBgain' = conversionFactor*DBgain + // DBped' = (conversionFactor*DBped*Dbgain � offset)/(conversionFactor*DBgain) + // = DBped - offset/DBgain' + // + + if (isLayer1) { + gain = gain * myVCalInfo.getConversionFactorL1(); + ped = ped - myVCalInfo.getOffsetL1() / gain; + } else { + gain = gain * myVCalInfo.getConversionFactor(); + ped = ped - myVCalInfo.getOffset() / gain; + } + + if (verbose_) + edm::LogInfo("SiPixelGainCalibScaler") << "post-change gain: " << gain << " pede:" << ped << std::endl; + + if constexpr (std::is_same_v) { + SiPixelGainCalibration_->setData(ped, gain, theSiPixelGainCalibration, false, false); + } else { + SiPixelGainCalibration_->setDataPedestal(ped, theSiPixelGainCalibration); + if ((row + 1) % numberOfRowsToAverageOver == 0) { // fill the column average after every ROC! + SiPixelGainCalibration_->setDataGain(gain, numberOfRowsToAverageOver, theSiPixelGainCalibration); + } + } + } // loop on rows + } // loop on columns + + typename PayloadType::Range outrange(theSiPixelGainCalibration.begin(), theSiPixelGainCalibration.end()); + if (!SiPixelGainCalibration_->put(d, outrange, ncols)) + edm::LogError("SiPixelGainCalibScaler") << "[SiPixelGainCalibScaler::analyze] detid already exists" << std::endl; + } // loop on DetIds + + // Write into DB + edm::LogInfo(" --- writing to DB!"); + edm::Service mydbservice; + if (!mydbservice.isAvailable()) { + edm::LogError("db service unavailable"); + return; + } else { + edm::LogInfo("DB service OK"); + } + + try { + if (mydbservice->isNewTagRequest(recordName_)) { + mydbservice->createNewIOV( + SiPixelGainCalibration_, mydbservice->beginOfTime(), mydbservice->endOfTime(), recordName_); + } else { + mydbservice->appendSinceTime(SiPixelGainCalibration_, mydbservice->currentTime(), recordName_); + } + edm::LogInfo(" --- all OK"); + } catch (const cond::Exception& er) { + edm::LogError("SiPixelGainCalibScaler") << er.what() << std::endl; + } catch (const std::exception& er) { + edm::LogError("SiPixelGainCalibScaler") << "caught std::exception " << er.what() << std::endl; + } catch (...) { + edm::LogError("SiPixelGainCalibScaler") << "Funny error" << std::endl; + } +} + +// ------------ method called once each job just before starting event loop ------------ +void SiPixelGainCalibScaler::beginJob() {} + +// ------------ method called once each job just after ending the event loop ------------ +void SiPixelGainCalibScaler::endJob() {} + +// ------------ method fills 'descriptions' with the allowed parameters for the module ------------ +void SiPixelGainCalibScaler::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + desc.add("record", "SiPixelGainCalibrationForHLTRcd"); + desc.add("isForHLT", true); + + edm::ParameterSetDescription vcalInfos; + vcalInfos.add("phase"); + vcalInfos.add("conversionFactor"); + vcalInfos.add("conversionFactorL1"); + vcalInfos.add("offset"); + vcalInfos.add("offsetL1"); + + std::vector tmp; + tmp.reserve(2); + { + edm::ParameterSet phase0VCal; + phase0VCal.addParameter("phase", 0); + phase0VCal.addParameter("conversionFactor", 65.); + phase0VCal.addParameter("conversionFactorL1", 65.); + phase0VCal.addParameter("offset", -414.); + phase0VCal.addParameter("offsetL1", -414.); + tmp.push_back(phase0VCal); + } + { + edm::ParameterSet phase1VCal; + phase1VCal.addParameter("phase", 1); + phase1VCal.addParameter("conversionFactor", 47.); + phase1VCal.addParameter("conversionFactorL1", 50.); + phase1VCal.addParameter("offset", -60.); + phase1VCal.addParameter("offsetL1", -670.); + tmp.push_back(phase1VCal); + } + desc.addVPSet("parameters", vcalInfos, tmp); + + desc.addUntracked("verbose", false); + + descriptions.add("siPixelGainCalibScaler", desc); +} + +//define this as a plug-in +DEFINE_FWK_MODULE(SiPixelGainCalibScaler); diff --git a/CondTools/SiPixel/test/SiPixelGainCalibScaler_cfg.py b/CondTools/SiPixel/test/SiPixelGainCalibScaler_cfg.py new file mode 100644 index 0000000000000..80eee74a06e02 --- /dev/null +++ b/CondTools/SiPixel/test/SiPixelGainCalibScaler_cfg.py @@ -0,0 +1,124 @@ +from __future__ import print_function +import FWCore.ParameterSet.Config as cms +import FWCore.ParameterSet.VarParsing as VarParsing + +process = cms.Process("Demo") + +## +## prepare options +## +options = VarParsing.VarParsing("analysis") + +options.register ('globalTag', + "103X_dataRun2_HLT_relval_v8",VarParsing.VarParsing.multiplicity.singleton, # singleton or list + VarParsing.VarParsing.varType.string, # string, int, or float + "GlobalTag") + +options.register ('forHLT', + True,VarParsing.VarParsing.multiplicity.singleton, # singleton or list + VarParsing.VarParsing.varType.bool, # string, int, or float + "is for SiPixelGainCalibrationForHLT") + +options.register ('firstRun', + 1,VarParsing.VarParsing.multiplicity.singleton, # singleton or list + VarParsing.VarParsing.varType.int, # string, int, or float + "first run to be processed") + +options.parseArguments() + +## +## MessageLogger +## +process.load('FWCore.MessageService.MessageLogger_cfi') +process.MessageLogger.categories.append("SiPixelGainCalibScaler") +process.MessageLogger.destinations = cms.untracked.vstring("cout") +process.MessageLogger.cout = cms.untracked.PSet( + threshold = cms.untracked.string("INFO"), + default = cms.untracked.PSet(limit = cms.untracked.int32(0)), + FwkReport = cms.untracked.PSet(limit = cms.untracked.int32(-1), + reportEvery = cms.untracked.int32(100000) + ), + SiPixelGainCalibScaler = cms.untracked.PSet( limit = cms.untracked.int32(-1)) + ) +process.MessageLogger.statistics.append('cout') + +process.load("Configuration.Geometry.GeometryRecoDB_cff") # Ideal geometry and interface +process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_cff') +from Configuration.AlCa.GlobalTag import GlobalTag +process.GlobalTag = GlobalTag(process.GlobalTag,options.globalTag, '') + +### Dirty trick to avoid geometry mismatches +process.trackerGeometryDB.applyAlignment = False + +## +## Selects the output record +## +MyRecord="" +if options.forHLT: + MyRecord="SiPixelGainCalibrationForHLTRcd" +else: + MyRecord="SiPixelGainCalibrationOfflineRcd" + +## +## Printing options +## +print("Using Global Tag:", process.GlobalTag.globaltag._value) +print("first run to be processed:",options.firstRun) +print("is for HLT? ","yes" if options.forHLT else "no!") +print("outputing on record: ",MyRecord) + +## +## Empty Source +## +process.maxEvents = cms.untracked.PSet( input = cms.untracked.int32(400000)) + +#################################################################### +# Empty source +#################################################################### +process.source = cms.Source("EmptySource", + firstRun = cms.untracked.uint32(options.firstRun), + numberEventsInRun = cms.untracked.uint32(1), # a number of events in single run + ) + +#################################################################### +# Analysis Module +#################################################################### +process.demo = cms.EDAnalyzer('SiPixelGainCalibScaler', + isForHLT = cms.bool(options.forHLT), + record = cms.string(MyRecord), + parameters = cms.VPSet( + cms.PSet( + conversionFactor = cms.double(65.), + conversionFactorL1 = cms.double(65.), + offset = cms.double(-414.), + offsetL1 = cms.double(-414.), + phase = cms.uint32(0) + ), + cms.PSet( + conversionFactor = cms.double(47.), + conversionFactorL1 = cms.double(50.), + offset = cms.double(-60.), + offsetL1 = cms.double(-670.), + phase = cms.uint32(1) + ) + ) + ) +## +## Database output service +## +process.load("CondCore.CondDB.CondDB_cfi") + +## +## Output database (in this case local sqlite file) +## +process.CondDB.connect = 'sqlite_file:TEST_modifiedGains_'+process.GlobalTag.globaltag._value+("_HLTGain" if options.forHLT else "_offlineGain")+".db" +process.PoolDBOutputService = cms.Service("PoolDBOutputService", + process.CondDB, + timetype = cms.untracked.string('runnumber'), + toPut = cms.VPSet(cms.PSet(record = cms.string(MyRecord), + tag = cms.string('scaledGains') + ) + ) + ) + +process.p = cms.Path(process.demo)