SiStrip/plugins/SiStripChannelGainFromDBMiscalibrator.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    CondTools/SiStrip
0004 // Class:      SiStripChannelGainFromDBMiscalibrator
0005 //
0006 /**\class SiStripChannelGainFromDBMiscalibrator SiStripChannelGainFromDBMiscalibrator.cc CondTools/SiStrip/plugins/SiStripChannelGainFromDBMiscalibrator.cc
0007
0008  Description: Class to miscalibrate a SiStrip Channel Gain payload from Database
0009
0010  Implementation:
0011      Read a SiStrip Channel Gain payload from DB (either central DB or sqlite file) and apply a miscalibration (either an offset / gaussian smearing or both)
0012      returns a local sqlite file with the same since of the original payload
0013 */
0014 //
0015 // Original Author:  Marco Musich
0016 //         Created:  Tue, 03 Oct 2017 12:57:34 GMT
0017 //
0018 //
0019
0020 // system include files
0021 #include <memory>
0022 #include <iostream>
0023
0024 // user include files
0025 #include "CLHEP/Random/RandGauss.h"
0026 #include "CalibFormats/SiStripObjects/interface/SiStripGain.h"
0027 #include "CalibTracker/Records/interface/SiStripGainRcd.h"
0028 #include "CommonTools/TrackerMap/interface/TrackerMap.h"
0029 #include "CondCore/DBOutputService/interface/PoolDBOutputService.h"
0030 #include "CondFormats/DataRecord/interface/SiStripApvGainRcd.h"
0031 #include "CondFormats/SiStripObjects/interface/SiStripApvGain.h"
0032 #include "CondFormats/SiStripObjects/interface/SiStripSummary.h"
0033 #include "CondTools/SiStrip/interface/SiStripMiscalibrateHelper.h"
0034 #include "DataFormats/SiStripDetId/interface/SiStripDetId.h"
0035 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0036 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0037 #include "FWCore/Framework/interface/ESHandle.h"
0038 #include "FWCore/Framework/interface/Event.h"
0039 #include "FWCore/Framework/interface/EventSetup.h"
0040 #include "FWCore/Framework/interface/Frameworkfwd.h"
0041 #include "FWCore/Framework/interface/MakerMacros.h"
0042 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0043 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0044 #include "FWCore/ServiceRegistry/interface/Service.h"
0045 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0046 //
0047 // class declaration
0048 //
0049
0050 class SiStripChannelGainFromDBMiscalibrator : public edm::one::EDAnalyzer<> {
0051 public:
0052   explicit SiStripChannelGainFromDBMiscalibrator(const edm::ParameterSet&);
0053   ~SiStripChannelGainFromDBMiscalibrator() override;
0054
0055   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0056
0057 private:
0058   void analyze(const edm::Event&, const edm::EventSetup&) override;
0059   SiStripApvGain getNewObject(const std::map<std::pair<uint32_t, int>, float>& theMap);
0060   void endJob() override;
0061
0062   // ----------member data ---------------------------
0063   const uint32_t m_printdebug;
0064   const std::string m_Record;
0065   const uint32_t m_gainType;
0066   const bool m_saveMaps;
0067   const std::vector<edm::ParameterSet> m_parameters;
0068   const edm::ESGetToken<SiStripGain, SiStripGainRcd> gainToken_;
0069   const edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> tTopoToken_;
0070
0071   std::unique_ptr<TrackerMap> scale_map;
0072   std::unique_ptr<TrackerMap> smear_map;
0073   std::unique_ptr<TrackerMap> ratio_map;
0074   std::unique_ptr<TrackerMap> old_payload_map;
0075   std::unique_ptr<TrackerMap> new_payload_map;
0076 };
0077
0078 //
0079 // constructors and destructor
0080 //
0081 SiStripChannelGainFromDBMiscalibrator::SiStripChannelGainFromDBMiscalibrator(const edm::ParameterSet& iConfig)
0082     : m_printdebug{iConfig.getUntrackedParameter<uint32_t>("printDebug", 1)},
0083       m_Record{iConfig.getUntrackedParameter<std::string>("record", "SiStripApvGainRcd")},
0084       m_gainType{iConfig.getUntrackedParameter<uint32_t>("gainType", 1)},
0085       m_saveMaps{iConfig.getUntrackedParameter<bool>("saveMaps", true)},
0086       m_parameters{iConfig.getParameter<std::vector<edm::ParameterSet> >("params")},
0087       gainToken_(esConsumes()),
0088       tTopoToken_(esConsumes()) {
0089   //now do what ever initialization is needed
0090
0091   std::string ss_gain = (m_gainType > 0) ? "G2" : "G1";
0092
0093   scale_map = std::make_unique<TrackerMap>("scale");
0094   scale_map->setTitle("Scale factor averaged by module");
0095   scale_map->setPalette(1);
0096
0097   smear_map = std::make_unique<TrackerMap>("smear");
0098   smear_map->setTitle("Smear factor averaged by module");
0099   smear_map->setPalette(1);
0100
0101   ratio_map = std::make_unique<TrackerMap>("ratio");
0102   ratio_map->setTitle("Average by module of the " + ss_gain + " Gain payload ratio (new/old)");
0103   ratio_map->setPalette(1);
0104
0105   new_payload_map = std::make_unique<TrackerMap>("new_payload");
0106   new_payload_map->setTitle("Tracker Map of Modified " + ss_gain + " Gain payload averaged by module");
0107   new_payload_map->setPalette(1);
0108
0109   old_payload_map = std::make_unique<TrackerMap>("old_payload");
0110   old_payload_map->setTitle("Tracker Map of Starting " + ss_gain + " Gain Payload averaged by module");
0111   old_payload_map->setPalette(1);
0112 }
0113
0114 SiStripChannelGainFromDBMiscalibrator::~SiStripChannelGainFromDBMiscalibrator() = default;
0115
0116 //
0117 // member functions
0118 //
0119
0120 // ------------ method called for each event  ------------
0121 void SiStripChannelGainFromDBMiscalibrator::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0122   using namespace edm;
0123
0124   const auto* const tTopo = &iSetup.getData(tTopoToken_);
0125
0126   std::vector<std::string> partitions;
0127
0128   // fill the list of partitions
0129   for (auto& thePSet : m_parameters) {
0130     const std::string partition(thePSet.getParameter<std::string>("partition"));
0131     // only if it is not yet in the list
0132     if (std::find(partitions.begin(), partitions.end(), partition) == partitions.end()) {
0133       partitions.push_back(partition);
0134     }
0135   }
0136
0137   std::map<sistripsummary::TrackerRegion, SiStripMiscalibrate::Smearings> mapOfSmearings;
0138
0139   for (auto& thePSet : m_parameters) {
0140     const std::string partition(thePSet.getParameter<std::string>("partition"));
0141     sistripsummary::TrackerRegion region = SiStripMiscalibrate::getRegionFromString(partition);
0142
0143     bool m_doScale(thePSet.getParameter<bool>("doScale"));
0144     bool m_doSmear(thePSet.getParameter<bool>("doSmear"));
0145     double m_scaleFactor(thePSet.getParameter<double>("scaleFactor"));
0146     double m_smearFactor(thePSet.getParameter<double>("smearFactor"));
0147
0148     SiStripMiscalibrate::Smearings params = SiStripMiscalibrate::Smearings();
0149     params.setSmearing(m_doScale, m_doSmear, m_scaleFactor, m_smearFactor);
0150     mapOfSmearings[region] = params;
0151   }
0152
0153   const auto& apvGain = iSetup.getData(gainToken_);
0154
0155   std::map<std::pair<uint32_t, int>, float> theMap, oldPayloadMap;
0156
0157   std::vector<uint32_t> detid;
0158   apvGain.getDetIds(detid);
0159   for (const auto& d : detid) {
0160     SiStripApvGain::Range range = apvGain.getRange(d, m_gainType);
0161     float nAPV = 0;
0162
0163     auto regions = SiStripMiscalibrate::getRegionsFromDetId(tTopo, d);
0164
0165     // sort by largest to smallest
0166     std::sort(regions.rbegin(), regions.rend());
0167
0168     SiStripMiscalibrate::Smearings params = SiStripMiscalibrate::Smearings();
0169
0170     for (unsigned int j = 0; j < regions.size(); j++) {
0171       bool checkRegion = (mapOfSmearings.count(regions[j]) != 0);
0172
0173       if (!checkRegion) {
0174         // if the subdetector is not in the list and there's no indication for the whole tracker, just use the default
0175         // i.e. no change
0176         continue;
0177       } else {
0178         params = mapOfSmearings[regions[j]];
0179         break;
0180       }
0181     }
0182
0183     scale_map->fill(d, params.m_scaleFactor);
0184     smear_map->fill(d, params.m_smearFactor);
0185
0186     for (int it = 0; it < range.second - range.first; it++) {
0187       nAPV += 1;
0188       float Gain = apvGain.getApvGain(it, range);
0189       std::pair<uint32_t, int> index = std::make_pair(d, nAPV);
0190
0191       oldPayloadMap[index] = Gain;
0192
0193       if (params.m_doScale) {
0194         Gain *= params.m_scaleFactor;
0195       }
0196
0197       if (params.m_doSmear) {
0198         float smearedGain = CLHEP::RandGauss::shoot(Gain, params.m_smearFactor);
0199         Gain = smearedGain;
0200       }
0201
0202       theMap[index] = Gain;
0203
0204     }  // loop over APVs
0205   }  // loop over DetIds
0206
0207   SiStripApvGain theAPVGains = this->getNewObject(theMap);
0208
0209   // make the payload ratio map
0210   uint32_t cachedId(0);
0211   SiStripMiscalibrate::Entry gain_ratio;
0212   SiStripMiscalibrate::Entry o_gain;
0213   SiStripMiscalibrate::Entry n_gain;
0214   unsigned int countDetIds(0);  // count DetIds to print
0215   for (const auto& element : theMap) {
0216     uint32_t DetId = element.first.first;
0217     int nAPV = element.first.second;
0218     float new_gain = element.second;
0219     float old_gain = oldPayloadMap[std::make_pair(DetId, nAPV)];
0220
0221     // flush the counters
0222     if (cachedId != 0 && DetId != cachedId) {
0223       ratio_map->fill(cachedId, gain_ratio.mean());
0224       old_payload_map->fill(cachedId, o_gain.mean());
0225       new_payload_map->fill(cachedId, n_gain.mean());
0226
0227       //auto test = new_payload_map.get()->smoduleMap;
0228
0229       gain_ratio.reset();
0230       o_gain.reset();
0231       n_gain.reset();
0232       countDetIds++;
0233     }
0234
0235     // printout for debug
0236     if (countDetIds < m_printdebug) {
0237       edm::LogPrint("SiStripChannelGainFromDBMiscalibrator")
0238           << "SiStripChannelGainFromDBMiscalibrator"
0239           << "::" << __FUNCTION__ << " detid " << DetId << " \t"
0240           << " APV " << nAPV << " \t new gain: " << new_gain << " \t old gain: " << old_gain << " \t" << std::endl;
0241     }
0242
0243     cachedId = DetId;
0244     gain_ratio.add(new_gain / old_gain);
0245     o_gain.add(old_gain);
0246     n_gain.add(new_gain);
0247   }
0248
0249   // write out the APVGains record
0250   edm::Service<cond::service::PoolDBOutputService> poolDbService;
0251
0252   if (poolDbService.isAvailable()) {
0253     if (poolDbService->isNewTagRequest(m_Record)) {
0254       poolDbService->createOneIOV(theAPVGains, poolDbService->currentTime(), m_Record);
0255     } else {
0256       poolDbService->appendOneIOV(theAPVGains, poolDbService->currentTime(), m_Record);
0257     }
0258   } else {
0259     throw std::runtime_error("PoolDBService required.");
0260   }
0261 }
0262
0263 // ------------ method called once each job just after ending the event loop  ------------
0264 void SiStripChannelGainFromDBMiscalibrator::endJob() {
0265   if (m_saveMaps) {
0266     std::string ss_gain = (m_gainType > 0) ? "G2" : "G1";
0267
0268     scale_map->save(true, 0, 0, ss_gain + "_gain_scale_map.pdf");
0269     scale_map->save(true, 0, 0, ss_gain + "_gain_scale_map.png");
0270
0271     smear_map->save(true, 0, 0, ss_gain + "_gain_smear_map.pdf");
0272     smear_map->save(true, 0, 0, ss_gain + "_gain_smear_map.png");
0273
0274     ratio_map->save(true, 0, 0, ss_gain + "_gain_ratio_map.pdf");
0275     ratio_map->save(true, 0, 0, ss_gain + "_gain_ratio_map.png");
0276
0277     auto range = SiStripMiscalibrate::getTruncatedRange(old_payload_map.get());
0278
0279     old_payload_map->save(true, range.first, range.second, "starting_" + ss_gain + "_gain_payload_map.pdf");
0280     old_payload_map->save(true, range.first, range.second, "starting_" + ss_gain + "_gain_payload_map.png");
0281
0282     range = SiStripMiscalibrate::getTruncatedRange(new_payload_map.get());
0283
0284     new_payload_map->save(true, range.first, range.second, "new_" + ss_gain + "_gain_payload_map.pdf");
0285     new_payload_map->save(true, range.first, range.second, "new_" + ss_gain + "_gain_payload_map.png");
0286   }
0287 }
0288
0289 //********************************************************************************//
0290 SiStripApvGain SiStripChannelGainFromDBMiscalibrator::getNewObject(
0291     const std::map<std::pair<uint32_t, int>, float>& theMap) {
0292   SiStripApvGain obj{};
0293
0294   std::vector<float> theSiStripVector;
0295   uint32_t PreviousDetId = 0;
0296   for (const auto& element : theMap) {
0297     uint32_t DetId = element.first.first;
0298     if (DetId != PreviousDetId) {
0299       if (!theSiStripVector.empty()) {
0300         SiStripApvGain::Range range(theSiStripVector.begin(), theSiStripVector.end());
0301         if (!obj.put(PreviousDetId, range))
0302           printf("Bug to put detId = %i\n", PreviousDetId);
0303       }
0304       theSiStripVector.clear();
0305       PreviousDetId = DetId;
0306     }
0307     theSiStripVector.push_back(element.second);
0308
0309     edm::LogInfo("SiStripChannelGainFromDBMiscalibrator")
0310         << " DetId: " << DetId << " APV:   " << element.first.second << " Gain:  " << element.second << std::endl;
0311   }
0312
0313   if (!theSiStripVector.empty()) {
0314     SiStripApvGain::Range range(theSiStripVector.begin(), theSiStripVector.end());
0315     if (!obj.put(PreviousDetId, range))
0316       printf("Bug to put detId = %i\n", PreviousDetId);
0317   }
0318
0319   return obj;
0320 }
0321
0322 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0323 void SiStripChannelGainFromDBMiscalibrator::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0324   edm::ParameterSetDescription desc;
0325
0326   desc.setComment(
0327       "Creates rescaled / smeared SiStrip Gain payload. Can be used for both G1 and G2."
0328       "PoolDBOutputService must be set up for 'SiStripApvGainRcd'.");
0329
0330   edm::ParameterSetDescription descScaler;
0331   descScaler.setComment(
0332       "ParameterSet specifying the Strip tracker partition to be scaled / smeared "
0333       "by a given factor.");
0334
0335   descScaler.add<std::string>("partition", "Tracker");
0336   descScaler.add<bool>("doScale", true);
0337   descScaler.add<bool>("doSmear", true);
0338   descScaler.add<double>("scaleFactor", 1.0);
0339   descScaler.add<double>("smearFactor", 1.0);
0340   desc.addVPSet("params", descScaler, std::vector<edm::ParameterSet>(1));
0341
0342   desc.addUntracked<unsigned int>("printDebug", 1);
0343   desc.addUntracked<std::string>("record", "SiStripApvGainRcd");
0344   desc.addUntracked<unsigned int>("gainType", 1);
0345   desc.addUntracked<bool>("saveMaps", true);
0346
0347   descriptions.add("scaleAndSmearSiStripGains", desc);
0348 }
0349
0350 /*--------------------------------------------------------------------*/
0351
0352 //define this as a plug-in
0353 DEFINE_FWK_MODULE(SiStripChannelGainFromDBMiscalibrator);