SiPhase2Tracker/plugins/SiPhase2BadStripChannelBuilder.cc

0001 // system include files
0002 #include <vector>
0003 #include <memory>
0004 #include <iostream>
0005 #include <fstream>
0006
0007 // user include files
0008 #include "CommonTools/ConditionDBWriter/interface/ConditionDBWriter.h"
0009 #include "CondFormats/SiStripObjects/interface/SiStripBadStrip.h"
0010 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0011 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0012 #include "DataFormats/Phase2TrackerDigi/interface/Phase2TrackerDigi.h"
0013 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0014 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0015 #include "FWCore/ParameterSet/interface/FileInPath.h"
0016 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0017 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0018 #include "FWCore/ServiceRegistry/interface/Service.h"
0019 #include "FWCore/Utilities/interface/Exception.h"
0020 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
0021 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0022 #include "Geometry/CommonTopologies/interface/PixelTopology.h"
0023 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0024 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0025 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0026
0027 #include <random>
0028
0029 /**
0030  * enum to decide which algorithm use to populate the conditions
0031  */
0032 namespace {
0033   enum badChannelAlgo { NAIVE = 1, RANDOM = 2, NONE = 99 };
0034 }
0035
0036 using Phase2TrackerGeomDetUnit = PixelGeomDetUnit;
0037
0038 class SiPhase2BadStripChannelBuilder : public ConditionDBWriter<SiStripBadStrip> {
0039 public:
0040   explicit SiPhase2BadStripChannelBuilder(const edm::ParameterSet&);
0041   ~SiPhase2BadStripChannelBuilder() override = default;
0042
0043   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0044
0045 private:
0046   std::unique_ptr<SiStripBadStrip> getNewObject() override;
0047
0048   void algoBeginRun(const edm::Run& run, const edm::EventSetup& es) override {
0049     if (!tTopo_) {
0050       tTopo_ = std::make_unique<TrackerTopology>(es.getData(topoToken_));
0051       const TrackerGeometry* tkGeom_ = &es.getData(geomToken_);
0052
0053       edm::LogInfo("SiPhase2BadStripChannelBuilder")
0054           << " There are " << tkGeom_->detUnits().size() << " modules in this geometry.";
0055
0056       for (auto const& det_u : tkGeom_->detUnits()) {
0057         const DetId detid = det_u->geographicalId();
0058         uint32_t rawId = detid.rawId();
0059         int subid = detid.subdetId();
0060         if (detid.det() == DetId::Detector::Tracker) {
0061           const Phase2TrackerGeomDetUnit* pixdet = dynamic_cast<const Phase2TrackerGeomDetUnit*>(det_u);
0062           assert(pixdet);
0063           LogDebug("SiPhase2BadStripChannelBuilder") << rawId << " is a " << subid << " det";
0064           if (subid == StripSubdetector::TOB || subid == StripSubdetector::TID) {
0065             if (tkGeom_->getDetectorType(rawId) == TrackerGeometry::ModuleType::Ph2PSS ||
0066                 tkGeom_->getDetectorType(rawId) == TrackerGeometry::ModuleType::Ph2SS) {
0067               theOTDets.push_back(pixdet);
0068             }
0069           }  // if it's a Strip module
0070         }  // if it's OT
0071       }  // if it's Tracker
0072     }  // loop of geomdets
0073   };
0074
0075   void algoAnalyze(const edm::Event& event, const edm::EventSetup& es) override {
0076     // deterministic seed from the event number
0077     // should not bias the result as the event number is already
0078     // assigned randomly-enough
0079     engine_.seed(event.id().event() + (event.id().luminosityBlock() << 10) + (event.id().run() << 20));
0080   }
0081
0082   std::map<unsigned short, unsigned short> clusterizeBadChannels(
0083       const std::vector<Phase2TrackerDigi::PackedDigiType>& maskedChannels);
0084
0085   // ----------member data ---------------------------
0086   std::unique_ptr<TrackerTopology> tTopo_;
0087   std::mt19937 engine_;
0088
0089   const edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> topoToken_;
0090   const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> geomToken_;
0091   const bool printdebug_;
0092   const unsigned int popConAlgo_;
0093   const float badComponentsFraction_;
0094   badChannelAlgo theBCAlgo_;
0095
0096   std::vector<const Phase2TrackerGeomDetUnit*> theOTDets;
0097 };
0098
0099 //__________________________________________________________________________________________________
0100 SiPhase2BadStripChannelBuilder::SiPhase2BadStripChannelBuilder(const edm::ParameterSet& iConfig)
0101     : ConditionDBWriter<SiStripBadStrip>(iConfig),
0102       topoToken_(esConsumes<edm::Transition::BeginRun>()),
0103       geomToken_(esConsumes<edm::Transition::BeginRun>()),
0104       printdebug_(iConfig.getUntrackedParameter<bool>("printDebug", false)),
0105       popConAlgo_(iConfig.getParameter<unsigned int>("popConAlgo")),
0106       badComponentsFraction_(iConfig.getParameter<double>("badComponentsFraction")) {
0107   if (badComponentsFraction_ > 1. || badComponentsFraction_ < 0.) {
0108     throw cms::Exception("Inconsistent configuration")
0109         << "[SiPhase2BadStripChannelBuilder::c'tor] the requested fraction of bad components is unphysical. \n";
0110   }
0111   theBCAlgo_ = static_cast<badChannelAlgo>(popConAlgo_);
0112 }
0113
0114 //__________________________________________________________________________________________________
0115 std::unique_ptr<SiStripBadStrip> SiPhase2BadStripChannelBuilder::getNewObject() {
0116   edm::LogInfo("SiPhase2BadStripChannelBuilder") << "... creating dummy SiStripBadStrip Data";
0117
0118   auto obj = std::make_unique<SiStripBadStrip>();
0119
0120   // early return with nullptr if fraction is ==0.
0121   if (badComponentsFraction_ == 0.f) {
0122     return obj;
0123   }
0124
0125   for (auto const& pixdet : theOTDets) {
0126     uint32_t rawId = pixdet->geographicalId().rawId();
0127     int subid = pixdet->geographicalId().subdetId();
0128
0129     const PixelTopology& topol(pixdet->specificTopology());
0130
0131     const int nrows = topol.nrows();
0132     const int ncols = topol.ncolumns();
0133
0134     LogDebug("SiPhase2BadStripChannelBuilder")
0135         << "DetId: " << rawId << " subdet: " << subid << " nrows: " << nrows << " ncols: " << ncols;
0136
0137     std::vector<unsigned int> theSiStripVector;
0138
0139     switch (theBCAlgo_) {
0140       case NAIVE: {
0141         LogDebug("SiPhase2BadStripChannelBuilder") << "using the NAIVE algorithm";
0142
0143         auto dis1 = std::uniform_int_distribution<>(0, nrows - 1);  // [0, nrows]
0144         auto dis2 = std::uniform_int_distribution<>(1, 10);         // [1, 10]
0145
0146         unsigned short firstBadStrip = std::floor(dis1(engine_));
0147         unsigned short NconsecutiveBadStrips = std::floor(dis2(engine_));
0148
0149         // if the interval exceeds the end of the module
0150         if (firstBadStrip + NconsecutiveBadStrips > nrows) {
0151           NconsecutiveBadStrips = nrows - firstBadStrip;
0152         }
0153
0154         unsigned int theBadStripRange;
0155         theBadStripRange = obj->encodePhase2(firstBadStrip, NconsecutiveBadStrips);
0156
0157         if (printdebug_)
0158           edm::LogInfo("SiPhase2BadStripChannelBuilder")
0159               << "detid " << rawId << " \t"
0160               << " firstBadStrip " << firstBadStrip << "\t "
0161               << " NconsecutiveBadStrips " << NconsecutiveBadStrips << "\t "
0162               << " packed integer " << std::hex << theBadStripRange << std::dec;
0163
0164         theSiStripVector.push_back(theBadStripRange);
0165         break;
0166       }
0167       case RANDOM: {
0168         LogDebug("SiPhase2BadStripChannelBuilder") << "using the RANDOM algorithm";
0169
0170         // auxilliary vector to check if the channels were already used
0171         std::vector<Phase2TrackerDigi::PackedDigiType> usedChannels;
0172
0173         size_t nmaxBadStrips = std::floor(nrows * ncols * badComponentsFraction_);
0174
0175         LogDebug("SiPhase2BadStripChannelBuilder")
0176             << __FUNCTION__ << " " << __LINE__ << " will mask: " << nmaxBadStrips << " strips";
0177
0178         auto disRows = std::uniform_int_distribution<>(0, nrows - 1);  // [0, nrows]
0179         auto disCols = std::uniform_int_distribution<>(0, ncols - 1);  // [0, ncols]
0180
0181         while (usedChannels.size() < nmaxBadStrips) {
0182           unsigned short badStripRow = std::floor(disRows(engine_));
0183           unsigned short badStripCol = std::floor(disCols(engine_));
0184
0185           const auto& badChannel = Phase2TrackerDigi::pixelToChannel(badStripRow, badStripCol);
0186
0187           LogDebug("SiPhase2BadStripChannelBuilder") << __FUNCTION__ << " " << __LINE__ << ": masking channel "
0188                                                      << badChannel << " (" << badStripRow << "," << badStripCol << ")";
0189
0190           if (std::find(usedChannels.begin(), usedChannels.end(), badChannel) == usedChannels.end()) {
0191             usedChannels.push_back(badChannel);
0192           }
0193         }
0194
0195         const auto badChannelsGroups = this->clusterizeBadChannels(usedChannels);
0196         // loop over the groups of bad strips
0197         for (const auto& [first, consec] : badChannelsGroups) {
0198           unsigned int theBadChannelsRange;
0199           theBadChannelsRange = obj->encodePhase2(first, consec);
0200
0201           if (printdebug_) {
0202             edm::LogInfo("SiPhase2BadStripChannelBuilder")
0203                 << "detid " << rawId << " \t"
0204                 << " firstBadStrip " << first << "\t "
0205                 << " NconsecutiveBadStrips " << consec << "\t "
0206                 << " packed integer " << std::hex << theBadChannelsRange << std::dec;
0207           }
0208           theSiStripVector.push_back(theBadChannelsRange);
0209         }
0210         break;
0211       }
0212       case NONE:
0213         [[fallthrough]];
0214       default:
0215         throw cms::Exception("Inconsistent configuration") << "Did not specifiy the right algorithm to be run";
0216     }
0217
0218     SiStripBadStrip::Range range(theSiStripVector.begin(), theSiStripVector.end());
0219     if (!obj->put(rawId, range))
0220       edm::LogError("SiPhase2BadStripChannelBuilder")
0221           << "[SiPhase2BadStripChannelBuilder::getNewObject] detid already exists";
0222
0223   }  // loop of geomdets
0224
0225   //End now write sistripbadChannel data in DB
0226   edm::Service<cond::service::PoolDBOutputService> mydbservice;
0227
0228   if (mydbservice.isAvailable()) {
0229     if (mydbservice->isNewTagRequest("SiStripBadStripRcd")) {
0230       mydbservice->createOneIOV<SiStripBadStrip>(*obj, mydbservice->beginOfTime(), "SiStripBadStripRcd");
0231     } else {
0232       mydbservice->appendOneIOV<SiStripBadStrip>(*obj, mydbservice->currentTime(), "SiStripBadStripRcd");
0233     }
0234   } else {
0235     edm::LogError("SiPhase2BadStripChannelBuilder") << "Service is unavailable";
0236   }
0237
0238   return obj;
0239 }
0240
0241 // poor-man clusterizing algorithm
0242 std::map<unsigned short, unsigned short> SiPhase2BadStripChannelBuilder::clusterizeBadChannels(
0243     const std::vector<Phase2TrackerDigi::PackedDigiType>& maskedChannels) {
0244   // Here we will store the result
0245   std::map<unsigned short, unsigned short> result{};
0246   std::map<int, std::string> printresult{};
0247
0248   // Sort and remove duplicates.
0249   std::set data(maskedChannels.begin(), maskedChannels.end());
0250
0251   // We will start the evaluation at the beginning of our data
0252   auto startOfSequence = data.begin();
0253
0254   // Find all sequences
0255   while (startOfSequence != data.end()) {
0256     // FInd first value that is not greate than one
0257     auto endOfSequence =
0258         std::adjacent_find(startOfSequence, data.end(), [](const auto& v1, const auto& v2) { return v2 != v1 + 1; });
0259     if (endOfSequence != data.end())
0260       std::advance(endOfSequence, 1);
0261
0262     auto consecutiveStrips = std::distance(startOfSequence, endOfSequence);
0263     result[*startOfSequence] = consecutiveStrips;
0264
0265     if (printdebug_) {
0266       // Build resulting string
0267       std::ostringstream oss{};
0268       bool writeDash = false;
0269       for (auto it = startOfSequence; it != endOfSequence; ++it) {
0270         oss << (writeDash ? "-" : "") << std::to_string(*it);
0271         writeDash = true;
0272       }
0273
0274       // Copy result to map
0275       for (auto it = startOfSequence; it != endOfSequence; ++it)
0276         printresult[*it] = oss.str();
0277     }
0278
0279     // Continue to search for the next sequence
0280     startOfSequence = endOfSequence;
0281   }
0282
0283   if (printdebug_) {
0284     // Show result on the screen. Or use the map in whichever way you want.
0285     for (const auto& [value, text] : printresult)
0286       edm::LogInfo("SiPhase2BadStripChannelBuilder") << std::left << std::setw(2) << value << " -> " << text << "\n";
0287   }
0288
0289   return result;
0290 }
0291
0292 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0293 void SiPhase2BadStripChannelBuilder::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0294   edm::ParameterSetDescription desc;
0295   desc.setComment("Module to build SiStripBadStrip Payloads for the Phase-2 Outer Tracker");
0296   ConditionDBWriter::fillPSetDescription(desc);  // inherited from mother class
0297   desc.addUntracked<bool>("printDebug", false)->setComment("maximum amount of print-outs");
0298   desc.add<unsigned int>("popConAlgo", 1)->setComment("algorithm to populate the payload: 1=NAIVE,2=RANDOM");
0299   desc.add<double>("badComponentsFraction", 0.01)->setComment("fraction of bad components to populate the payload");
0300   descriptions.addWithDefaultLabel(desc);
0301 }
0302
0303 #include "FWCore/PluginManager/interface/ModuleDef.h"
0304 #include "FWCore/Framework/interface/MakerMacros.h"
0305
0306 DEFINE_FWK_MODULE(SiPhase2BadStripChannelBuilder);