File indexing completed on 2023-03-17 11:27:50
0001
0002 #include <cmath>
0003 #include <fstream>
0004 #include <iostream>
0005 #include <map>
0006 #include <string>
0007 #include <vector>
0008
0009 #include "DataFormats/ForwardDetId/interface/ForwardSubdetector.h"
0010 #include "DataFormats/HGCRecHit/interface/HGCRecHitCollections.h"
0011 #include "DataFormats/DetId/interface/DetId.h"
0012
0013 #include "DQMServices/Core/interface/DQMEDAnalyzer.h"
0014 #include "DQMServices/Core/interface/DQMStore.h"
0015
0016 #include "FWCore/Framework/interface/Frameworkfwd.h"
0017 #include "FWCore/Framework/interface/ESHandle.h"
0018 #include "FWCore/Framework/interface/Event.h"
0019 #include "FWCore/Framework/interface/EventSetup.h"
0020 #include "FWCore/Utilities/interface/InputTag.h"
0021 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0022 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0023
0024 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0025 #include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
0026 #include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h"
0027 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0028
0029 #include "CLHEP/Units/GlobalSystemOfUnits.h"
0030 #include "TVector3.h"
0031
0032 class HGCalRecHitValidation : public DQMEDAnalyzer {
0033 public:
0034 struct energysum {
0035 energysum() { e15 = e25 = e50 = e100 = e250 = e1000 = 0.0; }
0036 double e15, e25, e50, e100, e250, e1000;
0037 };
0038
0039 struct HitsInfo {
0040 HitsInfo() {
0041 x = y = z = time = energy = phi = eta = 0.0;
0042 layer = 0;
0043 }
0044 float x, y, z, time, energy, phi, eta;
0045 float layer;
0046 };
0047
0048 explicit HGCalRecHitValidation(const edm::ParameterSet&);
0049 ~HGCalRecHitValidation() override = default;
0050
0051 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0052 void dqmBeginRun(const edm::Run&, const edm::EventSetup&) override;
0053 void bookHistograms(DQMStore::IBooker&, edm::Run const&, edm::EventSetup const&) override;
0054 void analyze(const edm::Event&, const edm::EventSetup&) override;
0055
0056 private:
0057 template <class T1, class T2>
0058 void recHitValidation(DetId& detId, int layer, const T1* geom, T2 it);
0059 void fillHitsInfo();
0060 void fillHitsInfo(HitsInfo& hits);
0061 void fillOccupancyMap(std::map<int, int>& OccupancyMap, int layer);
0062
0063
0064 const std::string nameDetector_;
0065 const int verbosity_;
0066 const edm::ESGetToken<HGCalGeometry, IdealGeometryRecord> geometry_token_;
0067 const edm::ESGetToken<HGCalGeometry, IdealGeometryRecord> geometry_beginRun_token_;
0068 const edm::EDGetTokenT<HGCRecHitCollection> recHitSource_;
0069 unsigned int layers_;
0070 int firstLayer_;
0071 std::map<int, int> OccupancyMap_plus;
0072 std::map<int, int> OccupancyMap_minus;
0073
0074 std::vector<MonitorElement*> EtaPhi_Plus_;
0075 std::vector<MonitorElement*> EtaPhi_Minus_;
0076 std::vector<MonitorElement*> energy_;
0077 std::vector<MonitorElement*> HitOccupancy_Plus_;
0078 std::vector<MonitorElement*> HitOccupancy_Minus_;
0079 MonitorElement* MeanHitOccupancy_Plus_;
0080 MonitorElement* MeanHitOccupancy_Minus_;
0081 };
0082
0083 HGCalRecHitValidation::HGCalRecHitValidation(const edm::ParameterSet& iConfig)
0084 : nameDetector_(iConfig.getParameter<std::string>("DetectorName")),
0085 verbosity_(iConfig.getUntrackedParameter<int>("Verbosity", 0)),
0086 geometry_token_(esConsumes(edm::ESInputTag("", nameDetector_))),
0087 geometry_beginRun_token_(esConsumes<HGCalGeometry, IdealGeometryRecord, edm::Transition::BeginRun>(
0088 edm::ESInputTag("", nameDetector_))),
0089 recHitSource_(consumes(iConfig.getParameter<edm::InputTag>("RecHitSource"))),
0090 firstLayer_(1) {}
0091
0092 void HGCalRecHitValidation::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0093 edm::ParameterSetDescription desc;
0094 desc.add<std::string>("DetectorName", "HGCalEESensitive");
0095 desc.add<edm::InputTag>("RecHitSource", edm::InputTag("HGCalRecHit", "HGCEERecHits"));
0096 desc.addUntracked<int>("Verbosity", 0);
0097 descriptions.add("hgcalRecHitValidationEE", desc);
0098 }
0099
0100 void HGCalRecHitValidation::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0101 OccupancyMap_plus.clear();
0102 OccupancyMap_minus.clear();
0103
0104 bool ok(true);
0105 unsigned int ntot(0), nused(0);
0106 const edm::ESHandle<HGCalGeometry>& geomHandle = iSetup.getHandle(geometry_token_);
0107 if (!geomHandle.isValid()) {
0108 edm::LogVerbatim("HGCalValidation") << "Cannot get valid HGCalGeometry "
0109 << "Object for " << nameDetector_;
0110 } else {
0111 const HGCalGeometry* geom0 = geomHandle.product();
0112 int geomType = ((geom0->topology().waferHexagon8()) ? 1 : ((geom0->topology().tileTrapezoid()) ? 2 : 0));
0113
0114 const edm::Handle<HGCRecHitCollection>& theRecHitContainers = iEvent.getHandle(recHitSource_);
0115 if (theRecHitContainers.isValid()) {
0116 if (verbosity_ > 0)
0117 edm::LogVerbatim("HGCalValidation")
0118 << nameDetector_ << " with " << theRecHitContainers->size() << " element(s)";
0119 for (const auto& it : *(theRecHitContainers.product())) {
0120 ntot++;
0121 nused++;
0122 DetId detId = it.id();
0123 int layer = ((geomType == 1) ? HGCSiliconDetId(detId).layer() : HGCScintillatorDetId(detId).layer());
0124 recHitValidation(detId, layer, geom0, &it);
0125 }
0126 } else {
0127 ok = false;
0128 edm::LogVerbatim("HGCalValidation") << "HGCRecHitCollection Handle "
0129 << "does not exist !!!";
0130 }
0131 }
0132 if (ok)
0133 fillHitsInfo();
0134 if (verbosity_ > 0)
0135 edm::LogVerbatim("HGCalValidation") << "Event " << iEvent.id().event() << " with " << ntot << " total and " << nused
0136 << " used recHits";
0137 }
0138
0139 template <class T1, class T2>
0140 void HGCalRecHitValidation::recHitValidation(DetId& detId, int layer, const T1* geom, T2 it) {
0141 const GlobalPoint& global = geom->getPosition(detId);
0142 double energy = it->energy();
0143
0144 float globalx = global.x();
0145 float globaly = global.y();
0146 float globalz = global.z();
0147
0148 HitsInfo hinfo;
0149 hinfo.energy = energy;
0150 hinfo.x = globalx;
0151 hinfo.y = globaly;
0152 hinfo.z = globalz;
0153 hinfo.layer = layer - firstLayer_;
0154 hinfo.phi = global.phi();
0155 hinfo.eta = global.eta();
0156
0157 if (verbosity_ > 1)
0158 edm::LogVerbatim("HGCalValidation") << "-------------------------- gx = " << globalx << " gy = " << globaly
0159 << " gz = " << globalz << " phi = " << hinfo.phi << " eta = " << hinfo.eta
0160 << " lay = " << hinfo.layer;
0161
0162 fillHitsInfo(hinfo);
0163
0164 if (hinfo.eta > 0)
0165 fillOccupancyMap(OccupancyMap_plus, hinfo.layer);
0166 else
0167 fillOccupancyMap(OccupancyMap_minus, hinfo.layer);
0168 }
0169
0170 void HGCalRecHitValidation::fillOccupancyMap(std::map<int, int>& OccupancyMap, int layer) {
0171 if (OccupancyMap.find(layer) != OccupancyMap.end())
0172 OccupancyMap[layer]++;
0173 else
0174 OccupancyMap[layer] = 1;
0175 }
0176
0177 void HGCalRecHitValidation::fillHitsInfo() {
0178 for (auto const& itr : OccupancyMap_plus) {
0179 int layer = itr.first;
0180 int occupancy = itr.second;
0181 HitOccupancy_Plus_.at(layer)->Fill(occupancy);
0182 }
0183
0184 for (auto const& itr : OccupancyMap_minus) {
0185 int layer = itr.first;
0186 int occupancy = itr.second;
0187 HitOccupancy_Minus_.at(layer)->Fill(occupancy);
0188 }
0189 }
0190
0191 void HGCalRecHitValidation::fillHitsInfo(HitsInfo& hits) {
0192 unsigned int ilayer = hits.layer;
0193 energy_.at(ilayer)->Fill(hits.energy);
0194
0195 EtaPhi_Plus_.at(ilayer)->Fill(hits.eta, hits.phi);
0196 EtaPhi_Minus_.at(ilayer)->Fill(hits.eta, hits.phi);
0197 }
0198
0199 void HGCalRecHitValidation::dqmBeginRun(const edm::Run&, const edm::EventSetup& iSetup) {
0200 const edm::ESHandle<HGCalGeometry>& geomHandle = iSetup.getHandle(geometry_beginRun_token_);
0201 if (!geomHandle.isValid()) {
0202 edm::LogVerbatim("HGCalValidation") << "Cannot get valid HGCalGeometry "
0203 << "Object for " << nameDetector_;
0204 } else {
0205 const HGCalGeometry* geom = geomHandle.product();
0206 const HGCalDDDConstants& hgcons_ = geom->topology().dddConstants();
0207 layers_ = hgcons_.layers(true);
0208 firstLayer_ = hgcons_.firstLayer();
0209 }
0210 }
0211
0212 void HGCalRecHitValidation::bookHistograms(DQMStore::IBooker& iB, edm::Run const&, edm::EventSetup const&) {
0213 iB.setCurrentFolder("HGCAL/HGCalRecHitsV/" + nameDetector_);
0214 std::ostringstream histoname;
0215 for (unsigned int il = 0; il < layers_; ++il) {
0216 int ilayer = firstLayer_ + static_cast<int>(il);
0217 auto istr1 = std::to_string(ilayer);
0218 while (istr1.size() < 2) {
0219 istr1.insert(0, "0");
0220 }
0221 histoname.str("");
0222 histoname << "HitOccupancy_Plus_layer_" << istr1;
0223 HitOccupancy_Plus_.push_back(iB.book1D(histoname.str().c_str(), "RecHitOccupancy_Plus", 100, 0, 10000));
0224 histoname.str("");
0225 histoname << "HitOccupancy_Minus_layer_" << istr1;
0226 HitOccupancy_Minus_.push_back(iB.book1D(histoname.str().c_str(), "RecHitOccupancy_Minus", 100, 0, 10000));
0227
0228 histoname.str("");
0229 histoname << "EtaPhi_Plus_"
0230 << "layer_" << istr1;
0231 EtaPhi_Plus_.push_back(iB.book2D(histoname.str().c_str(), "Occupancy", 31, 1.45, 3.0, 72, -CLHEP::pi, CLHEP::pi));
0232 histoname.str("");
0233 histoname << "EtaPhi_Minus_"
0234 << "layer_" << istr1;
0235 EtaPhi_Minus_.push_back(
0236 iB.book2D(histoname.str().c_str(), "Occupancy", 31, -3.0, -1.45, 72, -CLHEP::pi, CLHEP::pi));
0237
0238 histoname.str("");
0239 histoname << "energy_layer_" << istr1;
0240 energy_.push_back(iB.book1D(histoname.str().c_str(), "energy_", 500, 0, 1));
0241 }
0242
0243 histoname.str("");
0244 histoname << "SUMOfRecHitOccupancy_Plus";
0245 MeanHitOccupancy_Plus_ =
0246 iB.book1D(histoname.str().c_str(), "SUMOfRecHitOccupancy_Plus", layers_, -0.5, layers_ - 0.5);
0247 histoname.str("");
0248 histoname << "SUMOfRecHitOccupancy_Minus";
0249 MeanHitOccupancy_Minus_ =
0250 iB.book1D(histoname.str().c_str(), "SUMOfRecHitOccupancy_Minus", layers_, -0.5, layers_ - 0.5);
0251 }
0252
0253
0254
0255 #include "FWCore/Framework/interface/MakerMacros.h"
0256
0257
0258 DEFINE_FWK_MODULE(HGCalRecHitValidation);
