File indexing completed on 2024-05-10 02:21:34
0001
0002 #include <cmath>
0003 #include <iostream>
0004 #include <fstream>
0005 #include <vector>
0006 #include <map>
0007 #include <string>
0008
0009 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0010
0011 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0012
0013 #include "FWCore/Framework/interface/Frameworkfwd.h"
0014 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0015 #include "FWCore/Framework/interface/Event.h"
0016 #include "FWCore/Framework/interface/EventSetup.h"
0017 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0018 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0019 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0020 #include "FWCore/ServiceRegistry/interface/Service.h"
0021 #include "FWCore/Utilities/interface/InputTag.h"
0022
0023 #include "Geometry/HcalCommonData/interface/HcalHitRelabeller.h"
0024 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0025 #include "Geometry/HcalCommonData/interface/HcalDDDRecConstants.h"
0026 #include "Geometry/Records/interface/HcalRecNumberingRecord.h"
0027
0028 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0029 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0030
0031 #include <CLHEP/Geometry/Point3D.h>
0032 #include <CLHEP/Geometry/Transform3D.h>
0033 #include <CLHEP/Geometry/Vector3D.h>
0034 #include <CLHEP/Units/SystemOfUnits.h>
0035 #include <CLHEP/Units/PhysicalConstants.h>
0036
0037 #include "TH1D.h"
0038 #include "TH2D.h"
0039
0040 class HcalGeomCheck : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0041 public:
0042 struct hitsinfo {
0043 hitsinfo() { phi = eta = energy = time = 0.0; }
0044 double phi, eta, energy, time;
0045 };
0046
0047 explicit HcalGeomCheck(const edm::ParameterSet&);
0048 ~HcalGeomCheck() override = default;
0049 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0050
0051 protected:
0052 void beginJob() override;
0053 void beginRun(edm::Run const&, edm::EventSetup const&) override;
0054 void analyze(edm::Event const&, edm::EventSetup const&) override;
0055 void endRun(edm::Run const&, edm::EventSetup const&) override {}
0056
0057 private:
0058 void analyzeHits(int, const std::string&, const std::vector<PCaloHit>&);
0059
0060
0061 const std::string caloHitSource_;
0062 const int ietaMin_, ietaMax_, depthMax_;
0063 const double rmin_, rmax_, zmin_, zmax_;
0064 const int nbinR_, nbinZ_, verbosity_;
0065 const edm::EDGetTokenT<edm::PCaloHitContainer> tok_hits_;
0066 const edm::ESGetToken<HcalDDDRecConstants, HcalRecNumberingRecord> tok_HRNDC_;
0067 const HcalDDDRecConstants* hcons_;
0068
0069
0070 TH2D* h_RZ_;
0071 std::map<std::pair<int, int>, TH1D*> h_phi_;
0072 std::map<int, TH1D*> h_eta_;
0073 std::vector<TH1D*> h_E_, h_T_;
0074 };
0075
0076 HcalGeomCheck::HcalGeomCheck(const edm::ParameterSet& iConfig)
0077 : caloHitSource_(iConfig.getParameter<std::string>("caloHitSource")),
0078 ietaMin_(iConfig.getUntrackedParameter<int>("ietaMin", -41)),
0079 ietaMax_(iConfig.getUntrackedParameter<int>("ietaMax", 41)),
0080 depthMax_(iConfig.getUntrackedParameter<int>("depthMax", 7)),
0081 rmin_(iConfig.getUntrackedParameter<double>("rMin", 0.0)),
0082 rmax_(iConfig.getUntrackedParameter<double>("rMax", 5500.0)),
0083 zmin_(iConfig.getUntrackedParameter<double>("zMin", -12500.0)),
0084 zmax_(iConfig.getUntrackedParameter<double>("zMax", 12500.0)),
0085 nbinR_(iConfig.getUntrackedParameter<int>("nBinR", 550)),
0086 nbinZ_(iConfig.getUntrackedParameter<int>("nBinZ", 2500)),
0087 verbosity_(iConfig.getUntrackedParameter<int>("verbosity", 0)),
0088 tok_hits_(consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits", caloHitSource_))),
0089 tok_HRNDC_(esConsumes<HcalDDDRecConstants, HcalRecNumberingRecord, edm::Transition::BeginRun>()) {
0090 usesResource(TFileService::kSharedResource);
0091 }
0092
0093 void HcalGeomCheck::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0094 edm::ParameterSetDescription desc;
0095 desc.add<std::string>("caloHitSource", "HcalHits");
0096 desc.addUntracked<int>("ietaMin", -41);
0097 desc.addUntracked<int>("ietaMax", 41);
0098 desc.addUntracked<int>("depthMax", 7);
0099 desc.addUntracked<double>("rMin", 0.0);
0100 desc.addUntracked<double>("rMax", 5500.0);
0101 desc.addUntracked<double>("zMin", -12500.0);
0102 desc.addUntracked<double>("zMax", 12500.0);
0103 desc.addUntracked<int>("nBinR", 550);
0104 desc.addUntracked<int>("nBinZ", 250);
0105 desc.addUntracked<int>("verbosity", 0);
0106 descriptions.add("hcalGeomCheck", desc);
0107 }
0108
0109 void HcalGeomCheck::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0110
0111 const edm::Handle<edm::PCaloHitContainer>& theCaloHitContainer = iEvent.getHandle(tok_hits_);
0112 if (theCaloHitContainer.isValid()) {
0113 if (verbosity_ > 0)
0114 edm::LogVerbatim("HcalValidation") << " PcalohitItr = " << theCaloHitContainer->size();
0115
0116
0117 std::map<HcalDetId, hitsinfo> map_hits;
0118 for (auto const& hit : *(theCaloHitContainer.product())) {
0119 unsigned int id = hit.id();
0120 HcalDetId detId = HcalHitRelabeller::relabel(id, hcons_);
0121 double energy = hit.energy();
0122 double time = hit.time();
0123 int subdet = detId.subdet();
0124 int ieta = detId.ieta();
0125 int iphi = detId.iphi();
0126 int depth = detId.depth();
0127 std::pair<double, double> etaphi = hcons_->getEtaPhi(subdet, ieta, iphi);
0128 double rz = hcons_->getRZ(subdet, ieta, depth);
0129 if (verbosity_ > 2)
0130 edm::LogVerbatim("HcalValidation") << "i/p " << subdet << ":" << ieta << ":" << iphi << ":" << depth << " o/p "
0131 << etaphi.first << ":" << etaphi.second << ":" << rz;
0132 HepGeom::Point3D<float> gcoord = HepGeom::Point3D<float>(rz * cos(etaphi.second) / cosh(etaphi.first),
0133 rz * sin(etaphi.second) / cosh(etaphi.first),
0134 rz * tanh(etaphi.first));
0135 double tof = (gcoord.mag() * CLHEP::mm) / CLHEP::c_light;
0136 if (verbosity_ > 1)
0137 edm::LogVerbatim("HcalValidation")
0138 << "Detector " << subdet << " ieta = " << ieta << " iphi = " << iphi << " depth = " << depth
0139 << " positon = " << gcoord << " energy = " << energy << " time = " << time << ":" << tof;
0140 time -= tof;
0141 if (time < 0)
0142 time = 0;
0143 hitsinfo hinfo;
0144 if (map_hits.count(detId) != 0) {
0145 hinfo = map_hits[detId];
0146 } else {
0147 hinfo.phi = gcoord.getPhi();
0148 hinfo.eta = gcoord.getEta();
0149 hinfo.time = time;
0150 }
0151 hinfo.energy += energy;
0152 map_hits[detId] = hinfo;
0153
0154 h_RZ_->Fill(gcoord.z(), gcoord.rho());
0155 }
0156
0157
0158 for (auto const& hit : map_hits) {
0159 hitsinfo hinfo = hit.second;
0160 int subdet = hit.first.subdet();
0161 if (subdet > 0 && subdet < static_cast<int>(h_E_.size())) {
0162 int depth = hit.first.depth();
0163 int ieta = hit.first.ieta();
0164 int iphi = hit.first.iphi();
0165 if (verbosity_ > 1)
0166 edm::LogVerbatim("HGCalValidation")
0167 << " ---------------------- eta = " << ieta << ":" << hinfo.eta << " phi = " << iphi << ":" << hinfo.phi
0168 << " depth = " << depth << " E = " << hinfo.energy << " T = " << hinfo.time;
0169 h_E_[subdet]->Fill(hinfo.energy);
0170 h_T_[subdet]->Fill(hinfo.time);
0171 auto itr1 = h_phi_.find(std::pair<int, int>(ieta, depth));
0172 if (itr1 != h_phi_.end())
0173 itr1->second->Fill(iphi);
0174 auto itr2 = h_eta_.find(depth);
0175 if (itr2 != h_eta_.end())
0176 itr2->second->Fill(ieta);
0177 }
0178 }
0179 } else if (verbosity_ > 0) {
0180 edm::LogVerbatim("HcalValidation") << "PCaloHitContainer does not "
0181 << "exist for HCAL";
0182 }
0183 }
0184
0185
0186 void HcalGeomCheck::beginRun(const edm::Run&, const edm::EventSetup& iSetup) {
0187 const auto& pHRNDC = iSetup.getData(tok_HRNDC_);
0188 hcons_ = &pHRNDC;
0189 if (verbosity_ > 0)
0190 edm::LogVerbatim("HcalValidation") << "Obtain HcalDDDRecConstants from Event Setup";
0191 }
0192
0193 void HcalGeomCheck::beginJob() {
0194 if (verbosity_ > 2)
0195 edm::LogVerbatim("HcalValidation") << "HcalGeomCheck:: Enter beginJob";
0196 edm::Service<TFileService> fs;
0197 h_RZ_ = fs->make<TH2D>("RZ", "R vs Z", nbinZ_, zmin_, zmax_, nbinR_, rmin_, rmax_);
0198 if (verbosity_ > 2)
0199 edm::LogVerbatim("HcalValidation") << "HcalGeomCheck: booked scatterplot RZ";
0200 char name[20], title[100];
0201 for (int depth = 1; depth < depthMax_; ++depth) {
0202 for (int ieta = ietaMin_; ieta <= ietaMax_; ++ieta) {
0203 sprintf(name, "phi%d%d", ieta, depth);
0204 sprintf(title, "i#phi (i#eta = %d, depth = %d)", ieta, depth);
0205 h_phi_[std::pair<int, int>(ieta, depth)] = fs->make<TH1D>(name, title, 400, -20, 380);
0206 if (verbosity_ > 2)
0207 edm::LogVerbatim("HcalValidation") << "HcalGeomCheck: books " << title;
0208 }
0209 sprintf(name, "eta%d", depth);
0210 sprintf(title, "i#eta (depth = %d)", depth);
0211 h_eta_[depth] = fs->make<TH1D>(name, title, 100, -50, 50);
0212 if (verbosity_ > 2)
0213 edm::LogVerbatim("HcalValidation") << "HcalGeomCheck: books " << title;
0214 }
0215
0216 std::vector<std::string> dets = {"HCAL", "HB", "HE", "HF", "HO"};
0217 for (unsigned int ih = 0; ih < dets.size(); ++ih) {
0218 sprintf(name, "E_%s", dets[ih].c_str());
0219 sprintf(title, "Energy deposit in %s (MeV)", dets[ih].c_str());
0220 h_E_.emplace_back(fs->make<TH1D>(name, title, 1000, 0.0, 1.0));
0221 if (verbosity_ > 2)
0222 edm::LogVerbatim("HcalValidation") << "HcalGeomCheck: books " << title;
0223 sprintf(name, "T_%s", dets[ih].c_str());
0224 sprintf(title, "Time of hit in %s (ns)", dets[ih].c_str());
0225 h_T_.emplace_back(fs->make<TH1D>(name, title, 1000, 0.0, 200.0));
0226 if (verbosity_ > 2)
0227 edm::LogVerbatim("HcalValidation") << "HcalGeomCheck: books " << title;
0228 }
0229 }
0230
0231 #include "FWCore/Framework/interface/MakerMacros.h"
0232
0233 DEFINE_FWK_MODULE(HcalGeomCheck);