SiStripMonitorApproximateCluster/plugins/SiStripMonitorApproximateCluster.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    DQM/SiStripMonitorApproximateCluster
0004 // Class:      SiStripMonitorApproximateCluster
0005 //
0006 /**\class SiStripMonitorApproximateCluster SiStripMonitorApproximateCluster.cc DQM/SiStripMonitorApproximateCluster/plugins/SiStripMonitorApproximateCluster.cc
0007
0008  Description: Monitor SiStripApproximateClusters and on-demand compare properties with original SiStripClusters
0009
0010 */
0011 //
0012 // Original Author:  Marco Musich
0013 //         Created:  Thu, 08 Dec 2022 20:51:10 GMT
0014 //
0015 //
0016
0017 #include <string>
0018
0019 // user include files
0020 #include "DQMServices/Core/interface/DQMEDAnalyzer.h"
0021 #include "DQMServices/Core/interface/MonitorElement.h"
0022 #include "DataFormats/Common/interface/DetSet.h"
0023 #include "DataFormats/Common/interface/DetSetVectorNew.h"
0024 #include "DataFormats/SiStripCluster/interface/SiStripApproximateCluster.h"
0025 #include "DataFormats/SiStripCluster/interface/SiStripCluster.h"
0026 #include "FWCore/Framework/interface/Event.h"
0027 #include "FWCore/Framework/interface/Frameworkfwd.h"
0028 #include "FWCore/Framework/interface/MakerMacros.h"
0029 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0030 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0031 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h"
0032 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0033
0034 //
0035 // class declaration
0036 //
0037
0038 namespace siStripRawPrime {
0039   struct monitorApproxCluster {
0040   public:
0041     monitorApproxCluster()
0042         : h_barycenter_{nullptr}, h_width_{nullptr}, h_avgCharge_{nullptr}, h_isSaturated_{nullptr}, isBooked_{false} {}
0043
0044     void fill(const SiStripApproximateCluster& cluster) {
0045       h_barycenter_->Fill(cluster.barycenter());
0046       h_width_->Fill(cluster.width());
0047       h_avgCharge_->Fill(cluster.avgCharge());
0048       h_isSaturated_->Fill(cluster.isSaturated() ? 1 : -1);
0049     }
0050
0051     void book(dqm::implementation::DQMStore::IBooker& ibook, const std::string& folder) {
0052       ibook.setCurrentFolder(folder);
0053       h_barycenter_ =
0054           ibook.book1D("clusterBarycenter", "cluster barycenter;cluster barycenter;#clusters", 7680., 0., 7680.);
0055       h_width_ = ibook.book1D("clusterWidth", "cluster width;cluster width;#clusters", 128, -0.5, 127.5);
0056       h_avgCharge_ =
0057           ibook.book1D("clusterAvgCharge", "average strip charge;average strip charge;#clusters", 256, -0.5, 255.5);
0058       h_isSaturated_ = ibook.book1D("clusterSaturation", "cluster saturation;is saturated?;#clusters", 3, -1.5, 1.5);
0059       h_isSaturated_->getTH1F()->GetXaxis()->SetBinLabel(1, "Not saturated");
0060       h_isSaturated_->getTH1F()->GetXaxis()->SetBinLabel(3, "Saturated");
0061       isBooked_ = true;
0062     }
0063
0064     bool isBooked() { return isBooked_; }
0065
0066   private:
0067     dqm::reco::MonitorElement* h_barycenter_;
0068     dqm::reco::MonitorElement* h_width_;
0069     dqm::reco::MonitorElement* h_avgCharge_;
0070     dqm::reco::MonitorElement* h_isSaturated_;
0071     bool isBooked_;
0072   };
0073 }  // namespace siStripRawPrime
0074
0075 class SiStripMonitorApproximateCluster : public DQMEDAnalyzer {
0076 public:
0077   explicit SiStripMonitorApproximateCluster(const edm::ParameterSet&);
0078   ~SiStripMonitorApproximateCluster() override = default;
0079
0080   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0081
0082 private:
0083   void bookHistograms(DQMStore::IBooker&, edm::Run const&, edm::EventSetup const&) override;
0084
0085   void analyze(const edm::Event&, const edm::EventSetup&) override;
0086
0087   // ------------ member data ------------
0088   std::string folder_;
0089   bool compareClusters_;
0090   MonitorElement* h_nclusters_;
0091
0092   siStripRawPrime::monitorApproxCluster allClusters{};
0093   siStripRawPrime::monitorApproxCluster matchedClusters{};
0094   siStripRawPrime::monitorApproxCluster unMatchedClusters{};
0095
0096   // for comparisons
0097   MonitorElement* h_isMatched_{nullptr};
0098   MonitorElement* h_deltaBarycenter_{nullptr};
0099   MonitorElement* h_deltaSize_{nullptr};
0100   MonitorElement* h_deltaCharge_{nullptr};
0101   MonitorElement* h_deltaFirstStrip_{nullptr};
0102   MonitorElement* h_deltaEndStrip_{nullptr};
0103
0104   // Event Data
0105   edm::EDGetTokenT<edmNew::DetSetVector<SiStripApproximateCluster>> approxClustersToken_;
0106   edm::EDGetTokenT<edmNew::DetSetVector<SiStripCluster>> stripClustersToken_;
0107   const edmNew::DetSetVector<SiStripCluster>* stripClusterCollection_;
0108
0109   // Event Setup Data
0110   edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> tkGeomToken_;
0111 };
0112
0113 //
0114 // constructors and destructor
0115 //
0116 SiStripMonitorApproximateCluster::SiStripMonitorApproximateCluster(const edm::ParameterSet& iConfig)
0117     : folder_(iConfig.getParameter<std::string>("folder")),
0118       compareClusters_(iConfig.getParameter<bool>("compareClusters")),
0119       // Poducer name of input StripClusterCollection
0120       approxClustersToken_(consumes<edmNew::DetSetVector<SiStripApproximateCluster>>(
0121           iConfig.getParameter<edm::InputTag>("ApproxClustersProducer"))) {
0122   tkGeomToken_ = esConsumes();
0123   if (compareClusters_) {
0124     stripClustersToken_ =
0125         consumes<edmNew::DetSetVector<SiStripCluster>>(iConfig.getParameter<edm::InputTag>("ClustersProducer"));
0126   }
0127   stripClusterCollection_ = nullptr;
0128 }
0129
0130 //
0131 // member functions
0132 //
0133
0134 // ------------ method called for each event  ------------
0135 void SiStripMonitorApproximateCluster::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0136   using namespace edm;
0137
0138   const auto& tkGeom = &iSetup.getData(tkGeomToken_);
0139   const auto tkDets = tkGeom->dets();
0140
0141   // get collection of DetSetVector of clusters from Event
0142   edm::Handle<edmNew::DetSetVector<SiStripApproximateCluster>> approx_cluster_detsetvector;
0143   iEvent.getByToken(approxClustersToken_, approx_cluster_detsetvector);
0144   if (!approx_cluster_detsetvector.isValid()) {
0145     edm::LogError("SiStripMonitorApproximateCluster")
0146         << "SiStripApproximate cluster collection is not valid!" << std::endl;
0147
0148     // if approximate clusters collection not available, then early return
0149     return;
0150   }
0151
0152   // if requested to perform the comparison
0153   if (compareClusters_) {
0154     // get collection of DetSetVector of clusters from Event
0155     edm::Handle<edmNew::DetSetVector<SiStripCluster>> cluster_detsetvector;
0156     iEvent.getByToken(stripClustersToken_, cluster_detsetvector);
0157     if (!cluster_detsetvector.isValid()) {
0158       edm::LogError("SiStripMonitorApproximateCluster")
0159           << "Requested to perform comparison, but regular SiStrip cluster collection is not valid!" << std::endl;
0160       return;
0161     } else {
0162       stripClusterCollection_ = cluster_detsetvector.product();
0163     }
0164   }
0165
0166   int nApproxClusters{0};
0167   const edmNew::DetSetVector<SiStripApproximateCluster>* clusterCollection = approx_cluster_detsetvector.product();
0168
0169   for (const auto& detClusters : *clusterCollection) {
0170     edmNew::DetSet<SiStripCluster> strip_clusters_detset;
0171     const auto& detid = detClusters.detId();  // get the detid of the current detset
0172
0173     // starts here comaparison with regular clusters
0174     if (compareClusters_) {
0175       edmNew::DetSetVector<SiStripCluster>::const_iterator isearch =
0176           stripClusterCollection_->find(detid);  // search clusters of same detid
0177       strip_clusters_detset = (*isearch);
0178     }
0179
0180     for (const auto& cluster : detClusters) {
0181       nApproxClusters++;
0182
0183       // fill the full cluster collection
0184       allClusters.fill(cluster);
0185
0186       if (compareClusters_ && !strip_clusters_detset.empty()) {
0187         // build the converted cluster for the matching
0188         uint16_t nStrips{0};
0189         auto det = std::find_if(tkDets.begin(), tkDets.end(), [detid](auto& elem) -> bool {
0190           return (elem->geographicalId().rawId() == detid);
0191         });
0192         const StripTopology& p = dynamic_cast<const StripGeomDetUnit*>(*det)->specificTopology();
0193         nStrips = p.nstrips() - 1;
0194
0195         const auto convertedCluster = SiStripCluster(cluster, nStrips);
0196
0197         float distance{9999.};
0198         const SiStripCluster* closestCluster{nullptr};
0199         for (const auto& stripCluster : strip_clusters_detset) {
0200           // by construction the approximated cluster width has same
0201           // size as the original cluster
0202
0203           if (cluster.width() != stripCluster.size()) {
0204             continue;
0205           }
0206
0207           float deltaBarycenter = convertedCluster.barycenter() - stripCluster.barycenter();
0208           if (std::abs(deltaBarycenter) < distance) {
0209             closestCluster = &stripCluster;
0210             distance = std::abs(deltaBarycenter);
0211           }
0212         }
0213
0214         // Matching criteria:
0215         // - if exists a closest cluster in the DetId
0216         // - the size coincides with the original one
0217         if (closestCluster) {
0218           // comparisong plots
0219           h_deltaBarycenter_->Fill(closestCluster->barycenter() - convertedCluster.barycenter());
0220           h_deltaSize_->Fill(closestCluster->size() - convertedCluster.size());
0221           h_deltaCharge_->Fill(closestCluster->charge() - convertedCluster.charge());
0222           h_deltaFirstStrip_->Fill(closestCluster->firstStrip() - convertedCluster.firstStrip());
0223           h_deltaEndStrip_->Fill(closestCluster->endStrip() - convertedCluster.endStrip());
0224           // monitoring plots
0225           matchedClusters.fill(cluster);
0226           h_isMatched_->Fill(1);
0227         } else {
0228           // monitoring plots
0229           unMatchedClusters.fill(cluster);
0230           h_isMatched_->Fill(-1);
0231         }
0232       }  // if we're doing the comparison cluster by cluster
0233
0234     }  // loop on clusters in a detset
0235   }    // loop on the detset vector
0236   h_nclusters_->Fill(nApproxClusters);
0237 }
0238
0239 void SiStripMonitorApproximateCluster::bookHistograms(DQMStore::IBooker& ibook,
0240                                                       edm::Run const& run,
0241                                                       edm::EventSetup const& iSetup) {
0242   ibook.setCurrentFolder(folder_);
0243   h_nclusters_ = ibook.book1D("numberOfClusters", "total N. of clusters;N. of clusters;#clusters", 500., 0., 500000.);
0244   allClusters.book(ibook, folder_);
0245
0246   //  for comparisons
0247   if (compareClusters_) {
0248     // book monitoring for matche and unmatched clusters separately
0249     matchedClusters.book(ibook, fmt::format("{}/MatchedClusters", folder_));
0250     unMatchedClusters.book(ibook, fmt::format("{}/UnmatchedClusters", folder_));
0251
0252     ibook.setCurrentFolder(fmt::format("{}/ClusterComparisons", folder_));
0253     h_deltaBarycenter_ =
0254         ibook.book1D("deltaBarycenter", "#Delta barycenter;#Delta barycenter;cluster pairs", 101, -50.5, 50.5);
0255     h_deltaSize_ = ibook.book1D("deltaSize", "#Delta size;#Delta size;cluster pairs", 201, -100.5, 100.5);
0256     h_deltaCharge_ = ibook.book1D("deltaCharge", "#Delta charge;#Delta charge;cluster pairs", 401, -200.5, 200.5);
0257
0258     h_deltaFirstStrip_ =
0259         ibook.book1D("deltaFirstStrip", "#Delta FirstStrip; #Delta firstStrip;cluster pairs", 101, -50.5, 50.5);
0260     h_deltaEndStrip_ =
0261         ibook.book1D("deltaEndStrip", "#Delta EndStrip; #Delta endStrip; cluster pairs", 101, -50.5, 50.5);
0262
0263     h_isMatched_ = ibook.book1D("isClusterMatched", "cluster matching;is matched?;#clusters", 3, -1.5, 1.5);
0264     h_isMatched_->getTH1F()->GetXaxis()->SetBinLabel(1, "Not matched");
0265     h_isMatched_->getTH1F()->GetXaxis()->SetBinLabel(3, "Matched");
0266   }
0267 }
0268
0269 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0270 void SiStripMonitorApproximateCluster::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0271   edm::ParameterSetDescription desc;
0272   desc.setComment("Monitor SiStripApproximateCluster collection and compare with regular SiStrip clusters");
0273   desc.add<bool>("compareClusters", false)->setComment("if true, will compare with regualr Strip clusters");
0274   desc.add<edm::InputTag>("ApproxClustersProducer", edm::InputTag("hltSiStripClusters2ApproxClusters"))
0275       ->setComment("approxmate clusters collection");
0276   desc.add<edm::InputTag>("ClustersProducer", edm::InputTag("hltSiStripClusterizerForRawPrime"))
0277       ->setComment("regular clusters collection");
0278   desc.add<std::string>("folder", "SiStripApproximateClusters")->setComment("Top Level Folder");
0279   descriptions.addWithDefaultLabel(desc);
0280 }
0281
0282 // define this as a plug-in
0283 DEFINE_FWK_MODULE(SiStripMonitorApproximateCluster);