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
0004 #include "FWCore/Framework/interface/Event.h"
0005 #include "FWCore/Framework/interface/EventSetup.h"
0006 #include "FWCore/Framework/interface/MakerMacros.h"
0007
0008 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0009 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0010 #include "FWCore/Utilities/interface/InputTag.h"
0011 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0012 #include "FWCore/Utilities/interface/Exception.h"
0013
0014 #include "FWCore/ServiceRegistry/interface/Service.h"
0015 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0016
0017 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0018 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0019 #include "DataFormats/Math/interface/Point3D.h"
0020 #include "DataFormats/Math/interface/Vector3D.h"
0021
0022 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0023 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0024 #include "SimG4CMS/Calo/interface/CaloHitID.h"
0025
0026 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0027 #include "Geometry/HcalCommonData/interface/HcalDDDRecConstants.h"
0028 #include "Geometry/HcalCommonData/interface/HcalHitRelabeller.h"
0029 #include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"
0030 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0031 #include "Geometry/Records/interface/HcalRecNumberingRecord.h"
0032
0033 #include <TH2F.h>
0034 #include <cmath>
0035 #include <map>
0036 #include <memory>
0037 #include <string>
0038 #include <vector>
0039
0040 class HcalSimHitAnalysis : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0041 public:
0042 HcalSimHitAnalysis(const edm::ParameterSet& ps);
0043 ~HcalSimHitAnalysis() override = default;
0044 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0045
0046 protected:
0047 void beginJob() override;
0048 void analyze(edm::Event const&, edm::EventSetup const&) override;
0049 void beginRun(edm::Run const&, edm::EventSetup const&) override {}
0050 void endRun(edm::Run const&, edm::EventSetup const&) override {}
0051
0052 private:
0053 static const int ndets_ = 4;
0054 const std::string g4Label_, hcalHits_;
0055 const bool verbose_, testNumber_;
0056 const edm::EDGetTokenT<edm::PCaloHitContainer> tok_calo_;
0057 const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> tok_geom_;
0058 const edm::ESGetToken<HcalDDDRecConstants, HcalRecNumberingRecord> tok_HRNDC_;
0059 TH2F *poszp_[ndets_], *poszn_[ndets_];
0060 };
0061
0062 HcalSimHitAnalysis::HcalSimHitAnalysis(const edm::ParameterSet& ps)
0063 : g4Label_(ps.getUntrackedParameter<std::string>("moduleLabel", "g4SimHits")),
0064 hcalHits_(ps.getUntrackedParameter<std::string>("HitCollection", "HcalHits")),
0065 verbose_(ps.getUntrackedParameter<bool>("Verbose", false)),
0066 testNumber_(ps.getUntrackedParameter<bool>("TestNumber", true)),
0067 tok_calo_(consumes<edm::PCaloHitContainer>(edm::InputTag(g4Label_, hcalHits_))),
0068 tok_geom_(esConsumes<CaloGeometry, CaloGeometryRecord>()),
0069 tok_HRNDC_(esConsumes<HcalDDDRecConstants, HcalRecNumberingRecord>()) {
0070 usesResource(TFileService::kSharedResource);
0071
0072 edm::LogVerbatim("HitStudy") << "Module Label: " << g4Label_ << " Hits: " << hcalHits_ << " testNumber "
0073 << testNumber_;
0074 }
0075
0076 void HcalSimHitAnalysis::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0077 edm::ParameterSetDescription desc;
0078 desc.addUntracked<std::string>("ModuleLabel", "g4SimHits");
0079 desc.addUntracked<std::string>("HitCollection", "HcalHits");
0080 desc.addUntracked<bool>("Verbose", false);
0081 desc.addUntracked<bool>("TestNumber", true);
0082 descriptions.add("hcalSimHitAnalysis", desc);
0083 }
0084
0085 void HcalSimHitAnalysis::beginJob() {
0086 edm::Service<TFileService> tfile;
0087 if (!tfile.isAvailable())
0088 throw cms::Exception("BadConfig") << "TFileService unavailable: "
0089 << "please add it to config file";
0090 char name[20], title[120];
0091 std::string dets[ndets_] = {"HB", "HE", "HO", "HF"};
0092 int nx[ndets_] = {100, 100, 350, 160};
0093 double xlo[ndets_] = {0, -300, 0, -160};
0094 double xhi[ndets_] = {500, 300, 3500, 160};
0095 int ny[ndets_] = {100, 100, 50, 160};
0096 double ylo[ndets_] = {170, -300, 375, -160};
0097 double yhi[ndets_] = {370, 300, 425, 160};
0098 std::string xttl[ndets_] = {"|z| (cm)", "x (cm)", "|z| (cm)", "x (cm)"};
0099 std::string yttl[ndets_] = {"#rho (cm)", "y (cm)", "#rho (cm)", "y (cm)"};
0100 for (int i = 0; i < ndets_; i++) {
0101 sprintf(name, "poszp%d", i);
0102 sprintf(title, "%s+", dets[i].c_str());
0103 poszp_[i] = tfile->make<TH2F>(name, title, nx[i], xlo[i], xhi[i], ny[i], ylo[i], yhi[i]);
0104 poszp_[i]->GetXaxis()->SetTitle(xttl[i].c_str());
0105 poszp_[i]->GetYaxis()->SetTitle(yttl[i].c_str());
0106 sprintf(title, "%s-", dets[i].c_str());
0107 poszp_[i]->GetYaxis()->SetTitleOffset(1.2);
0108 poszp_[i]->Sumw2();
0109 sprintf(name, "poszn%d", i);
0110 poszn_[i] = tfile->make<TH2F>(name, title, nx[i], xlo[i], xhi[i], ny[i], ylo[i], yhi[i]);
0111 poszn_[i]->GetXaxis()->SetTitle(xttl[i].c_str());
0112 poszn_[i]->GetYaxis()->SetTitle(yttl[i].c_str());
0113 poszn_[i]->GetYaxis()->SetTitleOffset(1.2);
0114 poszn_[i]->Sumw2();
0115 }
0116 }
0117
0118 void HcalSimHitAnalysis::analyze(const edm::Event& e, const edm::EventSetup& iS) {
0119 if (verbose_)
0120 edm::LogVerbatim("HitStudy") << "Run = " << e.id().run() << " Event = " << e.id().event();
0121
0122
0123 const CaloGeometry* geo = &iS.getData(tok_geom_);
0124 const HcalGeometry* hgeom = static_cast<const HcalGeometry*>(geo->getSubdetectorGeometry(DetId::Hcal, HcalBarrel));
0125 const auto& pHRNDC = iS.getData(tok_HRNDC_);
0126 const HcalDDDRecConstants* hcons = &pHRNDC;
0127
0128 const edm::Handle<edm::PCaloHitContainer>& hitsCalo = e.getHandle(tok_calo_);
0129 bool getHits = (hitsCalo.isValid());
0130 uint32_t nhits = (getHits) ? hitsCalo->size() : 0;
0131 if (verbose_)
0132 edm::LogVerbatim("HitStudy") << "HcalSimHitAnalysis: Input flags Hits " << getHits << " with " << nhits << " hits";
0133 if (getHits) {
0134 std::vector<PCaloHit> hits;
0135 hits.insert(hits.end(), hitsCalo->begin(), hitsCalo->end());
0136 if (!hits.empty()) {
0137 std::map<HcalDetId, double> hitMap;
0138 for (auto hit : hits) {
0139 double edep = hit.energy();
0140 uint32_t id = hit.id();
0141 HcalDetId hid = (testNumber_) ? HcalHitRelabeller::relabel(id, hcons) : HcalDetId(id);
0142 auto it = hitMap.find(hid);
0143 if (it == hitMap.end()) {
0144 hitMap[hid] = edep;
0145 } else {
0146 (it->second) += edep;
0147 }
0148 }
0149
0150 for (auto it : hitMap) {
0151 HcalDetId id(it.first);
0152 GlobalPoint gpos = hgeom->getPosition(id);
0153 HcalSubdetector subdet = (id).subdet();
0154 int indx =
0155 ((subdet == HcalBarrel)
0156 ? 0
0157 : ((subdet == HcalEndcap) ? 1 : ((subdet == HcalOuter) ? 2 : ((subdet == HcalForward) ? 3 : -1))));
0158 if (verbose_)
0159 edm::LogVerbatim("HitStudy") << "HcalSimHitAnalysis: " << id << ":" << it.second << " at " << gpos
0160 << " subdet " << subdet << ":" << indx;
0161 if (indx >= 0) {
0162 double x = ((indx == 0) || (indx == 2)) ? std::abs(gpos.z()) : gpos.x();
0163 double y = ((indx == 0) || (indx == 2)) ? (gpos.perp()) : gpos.y();
0164 if (id.zside() >= 0)
0165 poszp_[indx]->Fill(x, y);
0166 else
0167 poszn_[indx]->Fill(x, y);
0168 }
0169 }
0170 }
0171 }
0172 }
0173
0174
0175 DEFINE_FWK_MODULE(HcalSimHitAnalysis);
