SiStrip/plugins/SiStripNoisesFromDBMiscalibrator.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    CondTools/SiStrip
0004 // Class:      SiStripNoisesFromDBMiscalibrator
0005 //
0006 /**\class SiStripNoisesFromDBMiscalibrator SiStripNoisesFromDBMiscalibrator.cc CondTools/SiStrip/plugins/SiStripNoisesFromDBMiscalibrator.cc
0007
0008  Description: Class to miscalibrate a SiStrip Noise payload from Database
0009
0010  Implementation:
0011      Read a SiStrip Noise 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 #include <fstream>
0024
0025 // user include files
0026 #include "CLHEP/Random/RandGauss.h"
0027 #include "CondTools/SiStrip/interface/SiStripMiscalibrateHelper.h"
0028 #include "CalibTracker/SiStripCommon/interface/SiStripDetInfoFileReader.h"
0029 #include "CommonTools/TrackerMap/interface/TrackerMap.h"
0030 #include "CondCore/DBOutputService/interface/PoolDBOutputService.h"
0031 #include "CondFormats/DataRecord/interface/SiStripNoisesRcd.h"
0032 #include "CondFormats/SiStripObjects/interface/SiStripNoises.h"
0033 #include "CondFormats/SiStripObjects/interface/SiStripSummary.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 //
0048 // class declaration
0049 //
0050
0051 class SiStripNoisesFromDBMiscalibrator : public edm::one::EDAnalyzer<> {
0052 public:
0053   explicit SiStripNoisesFromDBMiscalibrator(const edm::ParameterSet&);
0054   ~SiStripNoisesFromDBMiscalibrator() override = default;
0055
0056   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0057
0058 private:
0059   void analyze(const edm::Event&, const edm::EventSetup&) override;
0060   SiStripNoises getNewObject(const std::map<std::pair<uint32_t, int>, float>& theMap);
0061   SiStripNoises getNewObject_withDefaults(const std::map<std::pair<uint32_t, int>, float>& theMap,
0062                                           const float theDefault);
0063   void endJob() override;
0064
0065   // ----------member data ---------------------------
0066   const uint32_t m_printdebug;
0067   const bool m_perDetIDdebug;
0068   const bool m_fillDefaults;
0069   const bool m_saveMaps;
0070   const std::vector<edm::ParameterSet> m_parameters;
0071   const edm::FileInPath fp_;
0072   const edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> m_tTopoToken;
0073   const edm::ESGetToken<SiStripNoises, SiStripNoisesRcd> m_noiseToken;
0074
0075   std::unique_ptr<TrackerMap> scale_map;
0076   std::unique_ptr<TrackerMap> smear_map;
0077   std::unique_ptr<TrackerMap> ratio_map;
0078   std::unique_ptr<TrackerMap> old_payload_map;
0079   std::unique_ptr<TrackerMap> new_payload_map;
0080   std::unique_ptr<TrackerMap> missing_map;
0081 };
0082
0083 //
0084 // constructors and destructor
0085 //
0086 SiStripNoisesFromDBMiscalibrator::SiStripNoisesFromDBMiscalibrator(const edm::ParameterSet& iConfig)
0087     : m_printdebug{iConfig.getUntrackedParameter<uint32_t>("printDebug", 10)},
0088       m_perDetIDdebug{iConfig.getUntrackedParameter<bool>("perDetIDdebug", false)},
0089       m_fillDefaults{iConfig.getUntrackedParameter<bool>("fillDefaults", false)},
0090       m_saveMaps{iConfig.getUntrackedParameter<bool>("saveMaps", true)},
0091       m_parameters{iConfig.getParameter<std::vector<edm::ParameterSet> >("params")},
0092       fp_{iConfig.getUntrackedParameter<edm::FileInPath>("file",
0093                                                          edm::FileInPath(SiStripDetInfoFileReader::kDefaultFile))},
0094       m_tTopoToken(esConsumes()),
0095       m_noiseToken(esConsumes()) {
0096   //now do what ever initialization is needed
0097
0098   scale_map = std::make_unique<TrackerMap>("scale");
0099   scale_map->setTitle("Tracker Map of Scale factor averaged by module");
0100   scale_map->setPalette(1);
0101
0102   smear_map = std::make_unique<TrackerMap>("smear");
0103   smear_map->setTitle("Tracker Map of Smear factor averaged by module");
0104   smear_map->setPalette(1);
0105
0106   old_payload_map = std::make_unique<TrackerMap>("old_payload");
0107   old_payload_map->setTitle("Tracker Map of Starting Noise Payload averaged by module");
0108   old_payload_map->setPalette(1);
0109
0110   new_payload_map = std::make_unique<TrackerMap>("new_payload");
0111   new_payload_map->setTitle("Tracker Map of Modified Noise Payload averaged by module");
0112   new_payload_map->setPalette(1);
0113
0114   ratio_map = std::make_unique<TrackerMap>("ratio");
0115   ratio_map->setTitle("Tracker Map of Average by module of the payload ratio (new/old)");
0116   ratio_map->setPalette(1);
0117
0118   if (m_fillDefaults) {
0119     missing_map = std::make_unique<TrackerMap>("uncabled");
0120     missing_map->setTitle("Tracker Map of uncabled modules");
0121     missing_map->setPalette(1);
0122   }
0123 }
0124
0125 //
0126 // member functions
0127 //
0128
0129 // ------------ method called for each event  ------------
0130 void SiStripNoisesFromDBMiscalibrator::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0131   using namespace edm;
0132
0133   const auto tTopo = &iSetup.getData(m_tTopoToken);
0134
0135   std::vector<std::string> partitions;
0136
0137   // fill the list of partitions
0138   for (auto& thePSet : m_parameters) {
0139     const std::string partition(thePSet.getParameter<std::string>("partition"));
0140     // only if it is not yet in the list
0141     if (std::find(partitions.begin(), partitions.end(), partition) == partitions.end()) {
0142       partitions.push_back(partition);
0143     }
0144   }
0145
0146   std::map<sistripsummary::TrackerRegion, SiStripMiscalibrate::Smearings> mapOfSmearings;
0147
0148   for (auto& thePSet : m_parameters) {
0149     const std::string partition(thePSet.getParameter<std::string>("partition"));
0150     sistripsummary::TrackerRegion region = SiStripMiscalibrate::getRegionFromString(partition);
0151
0152     bool m_doScale(thePSet.getParameter<bool>("doScale"));
0153     bool m_doSmear(thePSet.getParameter<bool>("doSmear"));
0154     double m_scaleFactor(thePSet.getParameter<double>("scaleFactor"));
0155     double m_smearFactor(thePSet.getParameter<double>("smearFactor"));
0156
0157     SiStripMiscalibrate::Smearings params = SiStripMiscalibrate::Smearings();
0158     params.setSmearing(m_doScale, m_doSmear, m_scaleFactor, m_smearFactor);
0159     mapOfSmearings[region] = params;
0160   }
0161
0162   const auto& stripNoises = iSetup.getData(m_noiseToken);
0163
0164   std::map<std::pair<uint32_t, int>, float> theMap, oldPayloadMap;
0165
0166   std::vector<uint32_t> detid;
0167   stripNoises.getDetIds(detid);
0168   for (const auto& d : detid) {
0169     SiStripNoises::Range range = stripNoises.getRange(d);
0170
0171     auto regions = SiStripMiscalibrate::getRegionsFromDetId(tTopo, d);
0172
0173     // sort by largest to smallest
0174     std::sort(regions.rbegin(), regions.rend());
0175
0176     SiStripMiscalibrate::Smearings params = SiStripMiscalibrate::Smearings();
0177
0178     for (unsigned int j = 0; j < regions.size(); j++) {
0179       bool checkRegion = (mapOfSmearings.count(regions[j]) != 0);
0180
0181       if (!checkRegion) {
0182         // if the subdetector is not in the list and there's no indication for the whole tracker, just use the default
0183         // i.e. no change
0184         continue;
0185       } else {
0186         params = mapOfSmearings[regions[j]];
0187         break;
0188       }
0189     }
0190
0191     scale_map->fill(d, params.m_scaleFactor);
0192     smear_map->fill(d, params.m_smearFactor);
0193
0194     int nStrips = 0;
0195     for (int it = 0; it < (range.second - range.first) * 8 / 9; ++it) {
0196       auto noise = stripNoises.getNoise(it, range);
0197       std::pair<uint32_t, int> index = std::make_pair(d, nStrips);
0198
0199       oldPayloadMap[index] = noise;
0200
0201       if (params.m_doScale) {
0202         noise *= params.m_scaleFactor;
0203       }
0204
0205       if (params.m_doSmear) {
0206         float smearedNoise = CLHEP::RandGauss::shoot(noise, params.m_smearFactor);
0207         noise = smearedNoise;
0208       }
0209
0210       theMap[index] = noise;
0211
0212       nStrips += 1;
0213
0214     }  // loop over APVs
0215   }  // loop over DetIds
0216
0217   SiStripNoises theSiStripNoises{};
0218   if (!m_fillDefaults) {
0219     theSiStripNoises = this->getNewObject(theMap);
0220   } else {
0221     theSiStripNoises = this->getNewObject_withDefaults(theMap, -1.);
0222   }
0223
0224   // make the payload ratio map
0225   uint32_t cachedId(0);
0226   SiStripMiscalibrate::Entry noise_ratio;
0227   SiStripMiscalibrate::Entry o_noise;
0228   SiStripMiscalibrate::Entry n_noise;
0229   uint countDetIds(0);  // count DetIds to print
0230   uint countStrips(0);
0231   for (const auto& element : theMap) {
0232     countStrips++;
0233     uint32_t DetId = element.first.first;
0234     int nstrip = element.first.second;
0235     float new_noise = element.second;
0236     float old_noise = oldPayloadMap[std::make_pair(DetId, nstrip)];
0237
0238     // flush the counters
0239     if (cachedId != 0 && DetId != cachedId) {
0240       ratio_map->fill(cachedId, noise_ratio.mean());
0241       old_payload_map->fill(cachedId, o_noise.mean());
0242       new_payload_map->fill(cachedId, n_noise.mean());
0243
0244       //auto test = new_payload_map.get()->smoduleMap;
0245
0246       noise_ratio.reset();
0247       o_noise.reset();
0248       n_noise.reset();
0249       countDetIds++;
0250
0251       if (m_perDetIDdebug && (countDetIds < m_printdebug)) {
0252         edm::LogPrint("SiStripNoisesFromDBMiscalibrator")
0253             << "SiStripNoisesFromDBMiscalibrator"
0254             << "::" << __FUNCTION__ << " detid " << DetId << " \t"
0255             << " strip " << nstrip << " \t new <noise>: " << std::setw(5) << std::setprecision(2) << n_noise.mean()
0256             << " \t old <noise>: " << o_noise.mean() << " \t" << std::endl;
0257       }
0258     }
0259
0260     // printout for debug
0261     if ((countStrips < m_printdebug) && !m_perDetIDdebug) {
0262       edm::LogPrint("SiStripNoisesFromDBMiscalibrator")
0263           << "SiStripNoisesFromDBMiscalibrator"
0264           << "::" << __FUNCTION__ << " detid " << DetId << " \t"
0265           << " strip " << nstrip << " \t new noise: " << std::setw(5) << std::setprecision(2) << new_noise
0266           << " \t old noise: " << old_noise << " \t" << std::endl;
0267     }
0268
0269     cachedId = DetId;
0270     noise_ratio.add(new_noise / old_noise);
0271     o_noise.add(old_noise);
0272     n_noise.add(new_noise);
0273   }
0274
0275   // write out the SiStripNoises record
0276   edm::Service<cond::service::PoolDBOutputService> poolDbService;
0277
0278   if (poolDbService.isAvailable()) {
0279     if (poolDbService->isNewTagRequest("SiStripNoisesRcd")) {
0280       poolDbService->createOneIOV(theSiStripNoises, poolDbService->currentTime(), "SiStripNoisesRcd");
0281     } else {
0282       poolDbService->appendOneIOV(theSiStripNoises, poolDbService->currentTime(), "SiStripNoisesRcd");
0283     }
0284   } else {
0285     throw std::runtime_error("PoolDBService required.");
0286   }
0287 }
0288
0289 // ------------ method called once each job just after ending the event loop  ------------
0290 void SiStripNoisesFromDBMiscalibrator::endJob() {
0291   if (m_saveMaps) {
0292     scale_map->save(true, 0, 0, "noise_scale_map.pdf");
0293     scale_map->save(true, 0, 0, "noise_scale_map.png");
0294
0295     smear_map->save(true, 0, 0, "noise_smear_map.pdf");
0296     smear_map->save(true, 0, 0, "noise_smear_map.png");
0297
0298     ratio_map->save(true, 0, 0, "noise_ratio_map.pdf");
0299     ratio_map->save(true, 0, 0, "noise_ratio_map.png");
0300
0301     auto range = SiStripMiscalibrate::getTruncatedRange(old_payload_map.get());
0302
0303     old_payload_map->save(true, range.first, range.second, "starting_noise_payload_map.pdf");
0304     old_payload_map->save(true, range.first, range.second, "starting_noise_payload_map.png");
0305
0306     range = SiStripMiscalibrate::getTruncatedRange(new_payload_map.get());
0307
0308     new_payload_map->save(true, range.first, range.second, "new_noise_payload_map.pdf");
0309     new_payload_map->save(true, range.first, range.second, "new_noise_payload_map.png");
0310
0311     if (m_fillDefaults) {
0312       missing_map->save(true, 0, 0, "missing_map.pdf");
0313       missing_map->save(true, 0, 0, "missing_map.png");
0314     }
0315   }
0316 }
0317
0318 //********************************************************************************//
0319 SiStripNoises SiStripNoisesFromDBMiscalibrator::getNewObject_withDefaults(
0320     const std::map<std::pair<uint32_t, int>, float>& theMap, const float theDefault) {
0321   SiStripNoises obj{};
0322
0323   std::vector<uint32_t> missingDetIds;
0324
0325   const auto& reader = SiStripDetInfoFileReader::read(fp_.fullPath());
0326   const auto& DetInfos = reader.getAllData();
0327
0328   for (const auto& it : DetInfos) {
0329     const auto& nAPVs = it.second.nApvs;
0330     //Generate Noise for det detid
0331     bool isMissing(false);
0332     SiStripNoises::InputVector theSiStripVector;
0333     for (int t_strip = 0; t_strip < 128 * nAPVs; ++t_strip) {
0334       std::pair<uint32_t, int> index = std::make_pair(it.first, t_strip);
0335
0336       if (theMap.find(index) == theMap.end()) {
0337         LogDebug("SiStripNoisesFromDBMiscalibrator") << "detid " << it.first << " \t"
0338                                                      << " strip " << t_strip << " \t"
0339                                                      << " not found" << std::endl;
0340
0341         isMissing = true;
0342         obj.setData(theDefault, theSiStripVector);
0343
0344       } else {
0345         float noise = theMap.at(index);
0346         obj.setData(noise, theSiStripVector);
0347       }
0348     }
0349
0350     if (isMissing)
0351       missingDetIds.push_back(it.first);
0352
0353     if (!obj.put(it.first, theSiStripVector)) {
0354       edm::LogError("SiStripNoisesFromDBMiscalibrator")
0355           << "[SiStripNoisesFromDBMiscalibrator::analyze] detid already exists" << std::endl;
0356     }
0357   }
0358
0359   if (!missingDetIds.empty()) {
0360     // open output file
0361     std::stringstream name;
0362     name << "missing_modules.txt";
0363     std::ofstream* ofile = new std::ofstream(name.str(), std::ofstream::trunc);
0364     if (!ofile->is_open())
0365       throw "cannot open output file!";
0366     for (const auto& missing : missingDetIds) {
0367       edm::LogVerbatim("SiStripNoisesFromDBMiscalibrator") << missing << "  " << 1 << std::endl;
0368       (*ofile) << missing << " " << 1 << std::endl;
0369       missing_map->fill(missing, 1);
0370     }
0371
0372     ofile->close();
0373     delete ofile;
0374   }
0375
0376   return obj;
0377 }
0378
0379 //********************************************************************************//
0380 SiStripNoises SiStripNoisesFromDBMiscalibrator::getNewObject(const std::map<std::pair<uint32_t, int>, float>& theMap) {
0381   SiStripNoises obj{};
0382
0383   uint32_t PreviousDetId = 0;
0384   SiStripNoises::InputVector theSiStripVector;
0385   for (const auto& element : theMap) {
0386     uint32_t DetId = element.first.first;
0387     float noise = element.second;
0388
0389     if (DetId != PreviousDetId) {
0390       if (!theSiStripVector.empty()) {
0391         if (!obj.put(PreviousDetId, theSiStripVector)) {
0392           edm::LogError("SiStripNoisesFromDBMiscalibrator")
0393               << "[SiStripNoisesFromDBMiscalibrator::analyze] detid already exists" << std::endl;
0394         }
0395       }
0396
0397       theSiStripVector.clear();
0398       PreviousDetId = DetId;
0399     }
0400     obj.setData(noise, theSiStripVector);
0401   }
0402   return obj;
0403 }
0404
0405 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0406 void SiStripNoisesFromDBMiscalibrator::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0407   edm::ParameterSetDescription desc;
0408
0409   desc.setComment(
0410       "Creates rescaled / smeared SiStrip Noise payload."
0411       "PoolDBOutputService must be set up for 'SiSiStripNoisesRcd'.");
0412
0413   edm::ParameterSetDescription descScaler;
0414   descScaler.setComment(
0415       "ParameterSet specifying the Strip tracker partition to be scaled / smeared "
0416       "by a given factor.");
0417
0418   descScaler.add<std::string>("partition", "Tracker");
0419   descScaler.add<bool>("doScale", true);
0420   descScaler.add<bool>("doSmear", true);
0421   descScaler.add<double>("scaleFactor", 1.0);
0422   descScaler.add<double>("smearFactor", 1.0);
0423   desc.addVPSet("params", descScaler, std::vector<edm::ParameterSet>(1));
0424
0425   desc.addUntracked<unsigned int>("printDebug", 10);
0426   desc.addUntracked<bool>("perDetIDdebug", false);
0427   desc.addUntracked<bool>("fillDefaults", false);
0428   desc.addUntracked<bool>("saveMaps", true);
0429
0430   descriptions.add("scaleAndSmearSiStripNoises", desc);
0431 }
0432
0433 //define this as a plug-in
0434 DEFINE_FWK_MODULE(SiStripNoisesFromDBMiscalibrator);