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