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   for (const auto& it : SiStripDetInfoFileReader::read(fp_.fullPath()).getAllData()) {
0171     sprintf(hname, "h_%d", it.first);
0172     auto ref = ClusterPositionHistoMap.find(it.first);
0173     if (ref == ClusterPositionHistoMap.end()) {
0174       ClusterPositionHistoMap[it.first] =
0175           std::make_shared<TH1F>(hname, hname, it.second.nApvs * 128, -0.5, it.second.nApvs * 128 - 0.5);
0176     } else
0177       edm::LogError("SiStripQualityHotStripIdentifier")
0178           << " [SiStripQualityHotStripIdentifier::bookHistos] DetId " << it.first
0179           << " already found in map. Ignoring new data" << std::endl;
0180   }
0181 }
0182
0183 void SiStripQualityHotStripIdentifier::fillHisto(uint32_t detid, float value) {
0184   auto ref = ClusterPositionHistoMap.find(detid);
0185   if (ref != ClusterPositionHistoMap.end())
0186     ref->second->Fill(value);
0187   else
0188     edm::LogError("SiStripQualityHotStripIdentifier")
0189         << " [SiStripQualityHotStripIdentifier::fillHisto] Histogram not found in the list for DetId " << detid
0190         << " Ignoring data value " << value << std::endl;
0191 }
0192
0193 void SiStripQualityHotStripIdentifier::algoAnalyze(const edm::Event& e, const edm::EventSetup& eSetup) {
0194   edm::Handle<edm::DetSetVector<SiStripCluster> > dsv_SiStripCluster;
0195   e.getByLabel(Cluster_src_, dsv_SiStripCluster);
0196
0197   edm::Handle<reco::TrackCollection> trackCollection;
0198   if (tracksCollection_in_EventTree) {
0199     e.getByLabel(Track_src_, trackCollection);
0200     if (!trackCollection.isValid()) {
0201       edm::LogError("SiStripQualityHotStripIdentifier")
0202           << " [SiStripQualityHotStripIdentifier::algoAnalyze] missing trackCollection with label " << Track_src_
0203           << std::endl;
0204     }
0205   }
0206
0207   std::set<const void*> vPSiStripCluster;
0208   //Perform track study
0209   if (tracksCollection_in_EventTree) {
0210     int i = 0;
0211     for (const auto& track : *(trackCollection.product())) {
0212       LogTrace("SiStripQualityHotStripIdentifier")
0213           << "Track number " << i + 1 << "\n\tmomentum: " << track.momentum() << "\n\tPT: " << track.pt()
0214           << "\n\tvertex: " << track.vertex() << "\n\timpact parameter: " << track.d0()
0215           << "\n\tcharge: " << track.charge() << "\n\tnormalizedChi2: " << track.normalizedChi2() << "\n\tFrom EXTRA : "
0216           << "\n\t\touter PT " << track.outerPt() << std::endl;
0217
0218       //Loop on rechits
0219       for (auto const& recHit : track.recHits()) {
0220         if (!recHit->isValid()) {
0221           LogTrace("SiStripQualityHotStripIdentifier") << "\t\t Invalid Hit " << std::endl;
0222           continue;
0223         }
0224
0225         const SiStripRecHit2D* singleHit = dynamic_cast<const SiStripRecHit2D*>(recHit);
0226         const SiStripMatchedRecHit2D* matchedHit = dynamic_cast<const SiStripMatchedRecHit2D*>(recHit);
0227         const ProjectedSiStripRecHit2D* projectedHit = dynamic_cast<const ProjectedSiStripRecHit2D*>(recHit);
0228
0229         if (matchedHit) {
0230           vPSiStripCluster.insert((void*)&(matchedHit->monoCluster()));
0231           vPSiStripCluster.insert((void*)&(matchedHit->stereoCluster()));
0232         } else if (projectedHit) {
0233           vPSiStripCluster.insert((void*)&*(projectedHit->originalHit().cluster()));
0234         } else if (singleHit) {
0235           vPSiStripCluster.insert((void*)&*(singleHit->cluster()));
0236         } else {
0237           LogTrace("SiStripQualityHotStripIdentifier") << "NULL hit" << std::endl;
0238         }
0239       }
0240     }
0241   }
0242
0243   std::stringstream ss;
0244   //Loop on Det Clusters
0245   for (const auto& dSet : *dsv_SiStripCluster) {
0246     for (const auto& clus : dSet.data) {
0247       if (MinClusterWidth_ <= clus.amplitudes().size() && clus.amplitudes().size() <= MaxClusterWidth_) {
0248         if (std::find(vPSiStripCluster.begin(), vPSiStripCluster.end(), (void*)&clus) == vPSiStripCluster.end()) {
0249           if (edm::isDebugEnabled())
0250             ss << " adding cluster to histo for detid " << dSet.id << " with barycenter " << clus.barycenter()
0251                << std::endl;
0252           fillHisto(dSet.id, clus.barycenter());
0253         }
0254       }
0255     }
0256   }
0257   LogTrace("SiStripQualityHotStripIdentifier") << ss.str();
0258 }
0259
0260 #include "FWCore/Framework/interface/MakerMacros.h"
0261 DEFINE_FWK_MODULE(SiStripQualityHotStripIdentifier);