File indexing completed on 2024-04-06 12:29:50
0001 #include "FWCore/Framework/interface/Frameworkfwd.h"
0002 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0003 #include "FWCore/Framework/interface/Event.h"
0004 #include "FWCore/Framework/interface/EventSetup.h"
0005 #include "FWCore/Framework/interface/MakerMacros.h"
0006 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0007 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0008 #include "FWCore/Utilities/interface/InputTag.h"
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010 #include "FWCore/Utilities/interface/Exception.h"
0011
0012 #include "FWCore/ServiceRegistry/interface/Service.h"
0013 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0014
0015 #include "DataFormats/ForwardDetId/interface/HGCSiliconDetId.h"
0016 #include "DataFormats/ForwardDetId/interface/HGCScintillatorDetId.h"
0017
0018 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0019 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0020
0021 #include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h"
0022 #include "Geometry/CaloTopology/interface/HGCalTopology.h"
0023 #include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
0024 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0025
0026 #include <map>
0027 #include <string>
0028 #include <vector>
0029 #include <TH1D.h>
0030
0031 class HGcalHitCheck : public edm::one::EDAnalyzer<> {
0032 public:
0033 HGcalHitCheck(const edm::ParameterSet& ps);
0034 ~HGcalHitCheck() override {}
0035 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0036
0037 protected:
0038 void analyze(edm::Event const&, edm::EventSetup const&) override;
0039 void beginJob() override;
0040 void endJob() override {}
0041
0042 private:
0043 const std::string g4Label_, caloHitSource_, nameSense_, nameDetector_, tag_;
0044 const int layers_, verbosity_;
0045 const edm::EDGetTokenT<edm::PCaloHitContainer> tok_calo_;
0046 const edm::ESGetToken<HGCalGeometry, IdealGeometryRecord> geomToken_;
0047 bool histos_;
0048 TH1D *h_hits_, *h_hit1_, *h_hit2_;
0049 std::vector<TH1D*> h_hitL_, h_hitF_, h_hitP_;
0050 };
0051
0052 HGcalHitCheck::HGcalHitCheck(const edm::ParameterSet& ps)
0053 : g4Label_(ps.getParameter<std::string>("moduleLabel")),
0054 caloHitSource_(ps.getParameter<std::string>("caloHitSource")),
0055 nameSense_(ps.getParameter<std::string>("nameSense")),
0056 nameDetector_(ps.getParameter<std::string>("tag")),
0057 tag_(ps.getParameter<std::string>("nameDevice")),
0058 layers_(ps.getParameter<int>("layers")),
0059 verbosity_(ps.getParameter<int>("verbosity")),
0060 tok_calo_(consumes<edm::PCaloHitContainer>(edm::InputTag(g4Label_, caloHitSource_))),
0061 geomToken_(esConsumes<HGCalGeometry, IdealGeometryRecord>(edm::ESInputTag{"", nameSense_})),
0062 histos_(false) {
0063 edm::LogVerbatim("HitStudy") << "Test Hit ID for " << nameDetector_ << " using SimHits for " << nameSense_
0064 << " with module Label: " << g4Label_ << " Hits: " << caloHitSource_;
0065 }
0066
0067 void HGcalHitCheck::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0068 edm::ParameterSetDescription desc;
0069 desc.add<std::string>("moduleLabel", "g4SimHits");
0070 desc.add<std::string>("caloHitSource", "HGCHitsEE");
0071 desc.add<std::string>("nameSense", "HGCalEESensitive");
0072 desc.add<std::string>("nameDevice", "HGCal EE");
0073 desc.add<std::string>("tag", "DDD");
0074 desc.add<int>("layers", 26);
0075 desc.add<int>("verbosity", 0);
0076 descriptions.add("hgcalHitCheckEE", desc);
0077 }
0078
0079 void HGcalHitCheck::beginJob() {
0080 edm::Service<TFileService> fs;
0081 if (!fs.isAvailable()) {
0082 edm::LogVerbatim("HitStudy") << "TFileService unavailable: no histograms";
0083 } else {
0084 histos_ = true;
0085 char name[100], title[200];
0086 sprintf(name, "HitsL");
0087 sprintf(title, "Number of hits in %s for %s", nameSense_.c_str(), tag_.c_str());
0088 h_hits_ = fs->make<TH1D>(name, title, 1000, 0, 5000.);
0089 h_hits_->GetXaxis()->SetTitle(title);
0090 h_hits_->GetYaxis()->SetTitle("Hits");
0091 h_hits_->Sumw2();
0092 sprintf(name, "HitsF");
0093 sprintf(title, "Number of hits in %s for %s in Full Wafers or SiPM 2", nameSense_.c_str(), tag_.c_str());
0094 h_hit1_ = fs->make<TH1D>(name, title, 1000, 0, 5000.);
0095 h_hit1_->GetXaxis()->SetTitle(title);
0096 h_hit1_->GetYaxis()->SetTitle("Hits");
0097 h_hit1_->Sumw2();
0098 sprintf(name, "HitsP");
0099 sprintf(title, "Number of hits in %s for %s in Partial Wafers or SiPM 4", nameSense_.c_str(), tag_.c_str());
0100 h_hit2_ = fs->make<TH1D>(name, title, 1000, 0, 5000.);
0101 h_hit2_->GetXaxis()->SetTitle(title);
0102 h_hit2_->GetYaxis()->SetTitle("Hits");
0103 h_hit2_->Sumw2();
0104 for (int k = 0; k < layers_; ++k) {
0105 sprintf(name, "HitsL%d", k + 1);
0106 sprintf(title, "Number of hits in %s for %s in Layer %d", nameSense_.c_str(), tag_.c_str(), k + 1);
0107 h_hitL_.emplace_back(fs->make<TH1D>(name, title, 1000, 0, 5000.));
0108 h_hitL_.back()->GetXaxis()->SetTitle(title);
0109 h_hitL_.back()->GetYaxis()->SetTitle("Hits");
0110 h_hitL_.back()->Sumw2();
0111 sprintf(name, "HitsF%d", k + 1);
0112 sprintf(title,
0113 "Number of hits in %s for %s in Full Wafers or SiPM 2 of Layer %d",
0114 nameSense_.c_str(),
0115 tag_.c_str(),
0116 k + 1);
0117 h_hitF_.emplace_back(fs->make<TH1D>(name, title, 1000, 0, 5000.));
0118 h_hitF_.back()->GetXaxis()->SetTitle(title);
0119 h_hitF_.back()->GetYaxis()->SetTitle("Hits");
0120 h_hitF_.back()->Sumw2();
0121 sprintf(name, "HitsP%d", k + 1);
0122 sprintf(title,
0123 "Number of hits in %s for %s in Partial Wafers or SiPM 4 of Layer %d",
0124 nameSense_.c_str(),
0125 tag_.c_str(),
0126 k + 1);
0127 h_hitP_.emplace_back(fs->make<TH1D>(name, title, 1000, 0, 5000.));
0128 h_hitP_.back()->GetXaxis()->SetTitle(title);
0129 h_hitP_.back()->GetYaxis()->SetTitle("Hits");
0130 h_hitP_.back()->Sumw2();
0131 }
0132 }
0133 }
0134
0135 void HGcalHitCheck::analyze(const edm::Event& e, const edm::EventSetup& iS) {
0136 if (verbosity_ > 0)
0137 edm::LogVerbatim("HitStudy") << "Run = " << e.id().run() << " Event = " << e.id().event();
0138
0139
0140 const HGCalGeometry* geom = &iS.getData(geomToken_);
0141 const HGCalDDDConstants& hgc = geom->topology().dddConstants();
0142 const std::vector<DetId>& validIds = geom->getValidDetIds();
0143 edm::LogVerbatim("HitStudy") << "Detector " << nameSense_ << " with " << validIds.size() << " valid cells";
0144
0145 const edm::Handle<edm::PCaloHitContainer>& hitsCalo = e.getHandle(tok_calo_);
0146 bool getHits = (hitsCalo.isValid());
0147 uint32_t nhits = (getHits) ? hitsCalo->size() : 0;
0148 uint32_t wafer(0), tiles(0);
0149 if (verbosity_ > 1)
0150 edm::LogVerbatim("HitStudy") << "HGcalHitCheck: Input flags Hits " << getHits << " with " << nhits << " hits";
0151 if (histos_)
0152 h_hits_->Fill(nhits);
0153
0154 if (getHits) {
0155 std::vector<PCaloHit> hits;
0156 hits.insert(hits.end(), hitsCalo->begin(), hitsCalo->end());
0157 if (!hits.empty()) {
0158 for (auto hit : hits) {
0159 if (histos_) {
0160 if ((nameSense_ == "HGCalEESensitive") || (nameSense_ == "HGCalHESiliconSensitive")) {
0161 ++wafer;
0162 HGCSiliconDetId id(hit.id());
0163 int lay = id.layer();
0164 h_hitL_[lay - 1]->Fill(nhits);
0165 HGCalParameters::waferInfo info = hgc.waferInfo(lay, id.waferU(), id.waferV());
0166 if (info.part == HGCalTypes::WaferFull) {
0167 h_hit1_->Fill(nhits);
0168 h_hitF_[lay - 1]->Fill(nhits);
0169 } else {
0170 h_hit2_->Fill(nhits);
0171 h_hitP_[lay - 1]->Fill(nhits);
0172 }
0173 } else {
0174 ++tiles;
0175 HGCScintillatorDetId id(hit.id());
0176 int lay = id.layer();
0177 h_hitL_[lay - 1]->Fill(nhits);
0178 int sipm = id.sipm();
0179 if (sipm == 1) {
0180 h_hit2_->Fill(nhits);
0181 h_hitP_[lay - 1]->Fill(nhits);
0182 } else {
0183 h_hit1_->Fill(nhits);
0184 h_hitF_[lay - 1]->Fill(nhits);
0185 }
0186 }
0187 }
0188 }
0189 }
0190 }
0191 edm::LogVerbatim("HitStudy") << "Total hits = " << nhits << " Wafer DetIds = " << wafer << " Tile DetIds = " << tiles;
0192 }
0193
0194
0195 DEFINE_FWK_MODULE(HGcalHitCheck);
