SiStripQuality/plugins/SiStripQualityHotStripIdentifier.cc

0001 // system includes
0002 #include <fstream>
0003 #include <iostream>
0004 #include <memory>
0005 #include <sstream>
0006 #include <vector>
0007
0008 // user includes
0009 #include "CalibFormats/SiStripObjects/interface/SiStripQuality.h"
0010 #include "CalibTracker/Records/interface/SiStripQualityRcd.h"
0011 #include "CalibTracker/SiStripCommon/interface/SiStripDetInfoFileReader.h"
0012 #include "CalibTracker/SiStripQuality/interface/SiStripHotStripAlgorithmFromClusterOccupancy.h"
0013 #include "CalibTracker/SiStripQuality/interface/SiStripQualityHistos.h"
0014 #include "CommonTools/ConditionDBWriter/interface/ConditionDBWriter.h"
0015 #include "CondFormats/SiStripObjects/interface/SiStripBadStrip.h"
0016 #include "DataFormats/TrackReco/interface/Track.h"
0017 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0018 #include "DataFormats/TrackerRecHit2D/interface/ProjectedSiStripRecHit2D.h"
0019 #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h"
0020 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2D.h"
0021 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
0022 #include "FWCore/Framework/interface/ESWatcher.h"
0023 #include "FWCore/Framework/interface/EventSetup.h"
0024 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0025 #include "FWCore/ParameterSet/interface/FileInPath.h"
0026 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0027 #include "FWCore/Utilities/interface/Exception.h"
0028 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0029
0030 class TrackerTopology;
0031
0032 class SiStripQualityHotStripIdentifier : public ConditionDBWriter<SiStripBadStrip> {
0033 public:
0034   explicit SiStripQualityHotStripIdentifier(const edm::ParameterSet&);
0035   ~SiStripQualityHotStripIdentifier() override = default;
0036
0037 private:
0038   //Will be called at the beginning of the job
0039   void algoBeginJob(const edm::EventSetup&) override {}
0040   //Will be called at the beginning of each run in the job
0041   void algoBeginRun(const edm::Run&, const edm::EventSetup&) override;
0042   //Will be called at the beginning of each luminosity block in the run
0043   void algoBeginLuminosityBlock(const edm::LuminosityBlock&, const edm::EventSetup&) override { resetHistos(); }
0044   //Will be called at the end of the job
0045   void algoEndJob() override;
0046
0047   //Will be called at every event
0048   void algoAnalyze(const edm::Event&, const edm::EventSetup&) override;
0049
0050   std::unique_ptr<SiStripBadStrip> getNewObject() override;
0051
0052   void bookHistos();
0053   void resetHistos();
0054   void fillHisto(uint32_t detid, float value);
0055
0056 private:
0057   const std::string dataLabel_;
0058   const SiStripQuality* stripQuality_ = nullptr;
0059   const edm::ParameterSet conf_;
0060   const edm::FileInPath fp_;
0061   const edm::InputTag Cluster_src_;
0062   const edm::InputTag Track_src_;
0063   const bool tracksCollection_in_EventTree;
0064   const TrackerTopology* tTopo = nullptr;
0065
0066   unsigned short MinClusterWidth_, MaxClusterWidth_;
0067
0068   SiStrip::QualityHistosMap ClusterPositionHistoMap;
0069
0070   const edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> tTopoToken_;
0071   const edm::ESGetToken<SiStripQuality, SiStripQualityRcd> stripQualityToken_;
0072   edm::ESWatcher<SiStripQualityRcd> stripQualityWatcher_;
0073 };
0074
0075 SiStripQualityHotStripIdentifier::SiStripQualityHotStripIdentifier(const edm::ParameterSet& iConfig)
0076     : ConditionDBWriter<SiStripBadStrip>(iConfig),
0077       dataLabel_(iConfig.getUntrackedParameter<std::string>("dataLabel", "")),
0078       conf_(iConfig),
0079       fp_(iConfig.getUntrackedParameter<edm::FileInPath>("file",
0080                                                          edm::FileInPath(SiStripDetInfoFileReader::kDefaultFile))),
0081       Cluster_src_(iConfig.getParameter<edm::InputTag>("Cluster_src")),
0082       Track_src_(iConfig.getUntrackedParameter<edm::InputTag>("Track_src")),
0083       tracksCollection_in_EventTree(iConfig.getUntrackedParameter<bool>("RemoveTrackClusters", false)),
0084       tTopoToken_(esConsumes<edm::Transition::BeginRun>()),
0085       stripQualityToken_(esConsumes<edm::Transition::BeginRun>()) {
0086   edm::ParameterSet pset = iConfig.getUntrackedParameter<edm::ParameterSet>("ClusterSelection", edm::ParameterSet());
0087   MinClusterWidth_ = pset.getUntrackedParameter<uint32_t>("minWidth", 1);
0088   MaxClusterWidth_ = pset.getUntrackedParameter<uint32_t>("maxWidth", 1000);
0089
0090   bookHistos();
0091 }
0092
0093 std::unique_ptr<SiStripBadStrip> SiStripQualityHotStripIdentifier::getNewObject() {
0094   auto obj = std::make_unique<SiStripBadStrip>();
0095
0096   edm::ParameterSet parameters = conf_.getParameter<edm::ParameterSet>("AlgoParameters");
0097   std::string AlgoName = parameters.getParameter<std::string>("AlgoName");
0098   if (AlgoName == "SiStripHotStripAlgorithmFromClusterOccupancy") {
0099     edm::LogInfo("SiStripQualityHotStripIdentifier")
0100         << " [SiStripQualityHotStripIdentifier::getNewObject] call to SiStripHotStripAlgorithmFromClusterOccupancy"
0101         << std::endl;
0102
0103     SiStripHotStripAlgorithmFromClusterOccupancy theIdentifier(conf_, tTopo);
0104     theIdentifier.setProbabilityThreshold(parameters.getUntrackedParameter<double>("ProbabilityThreshold", 1.E-7));
0105     theIdentifier.setMinNumEntries(parameters.getUntrackedParameter<uint32_t>("MinNumEntries", 100));
0106     theIdentifier.setMinNumEntriesPerStrip(parameters.getUntrackedParameter<uint32_t>("MinNumEntriesPerStrip", 5));
0107
0108     const auto detInfo =
0109         SiStripDetInfoFileReader::read(edm::FileInPath{SiStripDetInfoFileReader::kDefaultFile}.fullPath());
0110     SiStripQuality* qobj = new SiStripQuality(detInfo);
0111     theIdentifier.extractBadStrips(qobj, ClusterPositionHistoMap, stripQuality_);
0112
0113     edm::LogInfo("SiStripQualityHotStripIdentifier")
0114         << " [SiStripQualityHotStripIdentifier::getNewObject] copy SiStripObject in SiStripBadStrip" << std::endl;
0115
0116     std::stringstream ss;
0117
0118     SiStripBadStrip::RegistryIterator rIter = qobj->getRegistryVectorBegin();
0119     SiStripBadStrip::RegistryIterator rIterEnd = qobj->getRegistryVectorEnd();
0120     for (; rIter != rIterEnd; ++rIter) {
0121       SiStripBadStrip::Range range(qobj->getDataVectorBegin() + rIter->ibegin,
0122                                    qobj->getDataVectorBegin() + rIter->iend);
0123       if (!obj->put(rIter->detid, range))
0124         edm::LogError("SiStripQualityHotStripIdentifier")
0125             << "[SiStripQualityHotStripIdentifier::getNewObject] detid already exists" << std::endl;
0126     }
0127     edm::LogInfo("SiStripQualityHotStripIdentifier")
0128         << " [SiStripQualityHotStripIdentifier::getNewObject] " << ss.str() << std::endl;
0129
0130   } else {
0131     edm::LogError("SiStripQualityHotStripIdentifier")
0132         << " [SiStripQualityHotStripIdentifier::getNewObject] call for a unknow HotStrip identification algoritm"
0133         << std::endl;
0134
0135     std::vector<uint32_t> a;
0136     SiStripBadStrip::Range range(a.begin(), a.end());
0137     if (!obj->put(0xFFFFFFFF, range))
0138       edm::LogError("SiStripQualityHotStripIdentifier")
0139           << "[SiStripQualityHotStripIdentifier::getNewObject] detid already exists" << std::endl;
0140   }
0141
0142   return obj;
0143 }
0144
0145 void SiStripQualityHotStripIdentifier::algoBeginRun(const edm::Run& run, const edm::EventSetup& iSetup) {
0146   tTopo = &iSetup.getData(tTopoToken_);
0147
0148   resetHistos();
0149
0150   if (stripQualityWatcher_.check(iSetup)) {
0151     stripQuality_ = &iSetup.getData(stripQualityToken_);
0152   }
0153 }
0154
0155 void SiStripQualityHotStripIdentifier::algoEndJob() {
0156   //Clear map
0157   ClusterPositionHistoMap.clear();
0158 }
0159
0160 void SiStripQualityHotStripIdentifier::resetHistos() {
0161   edm::LogInfo("SiStripQualityHotStripIdentifier") << " [SiStripQualityHotStripIdentifier::resetHistos] " << std::endl;
0162   for (const auto& it : ClusterPositionHistoMap) {
0163     it.second->Reset();
0164   }
0165 }
0166
0167 void SiStripQualityHotStripIdentifier::bookHistos() {
0168   edm::LogInfo("SiStripQualityHotStripIdentifier") << " [SiStripQualityHotStripIdentifier::bookHistos] " << std::endl;
0169   char hname[1024];
0170   const SiStripDetInfo& info = SiStripDetInfoFileReader::read(fp_.fullPath());
0171   for (const auto& it : info.getAllData()) {
0172     sprintf(hname, "h_%d", it.first);
0173     auto ref = ClusterPositionHistoMap.find(it.first);
0174     if (ref == ClusterPositionHistoMap.end()) {
0175       ClusterPositionHistoMap[it.first] =
0176           std::make_shared<TH1F>(hname, hname, it.second.nApvs * 128, -0.5, it.second.nApvs * 128 - 0.5);
0177     } else
0178       edm::LogError("SiStripQualityHotStripIdentifier")
0179           << " [SiStripQualityHotStripIdentifier::bookHistos] DetId " << it.first
0180           << " already found in map. Ignoring new data" << std::endl;
0181   }
0182 }
0183
0184 void SiStripQualityHotStripIdentifier::fillHisto(uint32_t detid, float value) {
0185   auto ref = ClusterPositionHistoMap.find(detid);
0186   if (ref != ClusterPositionHistoMap.end())
0187     ref->second->Fill(value);
0188   else
0189     edm::LogError("SiStripQualityHotStripIdentifier")
0190         << " [SiStripQualityHotStripIdentifier::fillHisto] Histogram not found in the list for DetId " << detid
0191         << " Ignoring data value " << value << std::endl;
0192 }
0193
0194 void SiStripQualityHotStripIdentifier::algoAnalyze(const edm::Event& e, const edm::EventSetup& eSetup) {
0195   edm::Handle<edm::DetSetVector<SiStripCluster> > dsv_SiStripCluster;
0196   e.getByLabel(Cluster_src_, dsv_SiStripCluster);
0197
0198   edm::Handle<reco::TrackCollection> trackCollection;
0199   if (tracksCollection_in_EventTree) {
0200     e.getByLabel(Track_src_, trackCollection);
0201     if (!trackCollection.isValid()) {
0202       edm::LogError("SiStripQualityHotStripIdentifier")
0203           << " [SiStripQualityHotStripIdentifier::algoAnalyze] missing trackCollection with label " << Track_src_
0204           << std::endl;
0205     }
0206   }
0207
0208   std::set<const void*> vPSiStripCluster;
0209   //Perform track study
0210   if (tracksCollection_in_EventTree) {
0211     int i = 0;
0212     for (const auto& track : *(trackCollection.product())) {
0213       LogTrace("SiStripQualityHotStripIdentifier")
0214           << "Track number " << i + 1 << "\n\tmomentum: " << track.momentum() << "\n\tPT: " << track.pt()
0215           << "\n\tvertex: " << track.vertex() << "\n\timpact parameter: " << track.d0()
0216           << "\n\tcharge: " << track.charge() << "\n\tnormalizedChi2: " << track.normalizedChi2() << "\n\tFrom EXTRA : "
0217           << "\n\t\touter PT " << track.outerPt() << std::endl;
0218
0219       //Loop on rechits
0220       for (auto const& recHit : track.recHits()) {
0221         if (!recHit->isValid()) {
0222           LogTrace("SiStripQualityHotStripIdentifier") << "\t\t Invalid Hit " << std::endl;
0223           continue;
0224         }
0225
0226         const SiStripRecHit2D* singleHit = dynamic_cast<const SiStripRecHit2D*>(recHit);
0227         const SiStripMatchedRecHit2D* matchedHit = dynamic_cast<const SiStripMatchedRecHit2D*>(recHit);
0228         const ProjectedSiStripRecHit2D* projectedHit = dynamic_cast<const ProjectedSiStripRecHit2D*>(recHit);
0229
0230         if (matchedHit) {
0231           vPSiStripCluster.insert((void*)&(matchedHit->monoCluster()));
0232           vPSiStripCluster.insert((void*)&(matchedHit->stereoCluster()));
0233         } else if (projectedHit) {
0234           vPSiStripCluster.insert((void*)&*(projectedHit->originalHit().cluster()));
0235         } else if (singleHit) {
0236           vPSiStripCluster.insert((void*)&*(singleHit->cluster()));
0237         } else {
0238           LogTrace("SiStripQualityHotStripIdentifier") << "NULL hit" << std::endl;
0239         }
0240       }
0241     }
0242   }
0243
0244   std::stringstream ss;
0245   //Loop on Det Clusters
0246   for (const auto& dSet : *dsv_SiStripCluster) {
0247     for (const auto& clus : dSet.data) {
0248       if (MinClusterWidth_ <= clus.amplitudes().size() && clus.amplitudes().size() <= MaxClusterWidth_) {
0249         if (std::find(vPSiStripCluster.begin(), vPSiStripCluster.end(), (void*)&clus) == vPSiStripCluster.end()) {
0250           if (edm::isDebugEnabled())
0251             ss << " adding cluster to histo for detid " << dSet.id << " with barycenter " << clus.barycenter()
0252                << std::endl;
0253           fillHisto(dSet.id, clus.barycenter());
0254         }
0255       }
0256     }
0257   }
0258   LogTrace("SiStripQualityHotStripIdentifier") << ss.str();
0259 }
0260
0261 #include "FWCore/Framework/interface/MakerMacros.h"
0262 DEFINE_FWK_MODULE(SiStripQualityHotStripIdentifier);