File indexing completed on 2023-03-17 11:27:50
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/ForwardDetId/interface/HFNoseDetId.h"
0012 #include "DataFormats/ForwardDetId/interface/HGCSiliconDetId.h"
0013 #include "DataFormats/ForwardDetId/interface/HGCScintillatorDetId.h"
0014
0015 #include "FWCore/Framework/interface/Frameworkfwd.h"
0016 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0017 #include "FWCore/Framework/interface/Event.h"
0018 #include "FWCore/Framework/interface/EventSetup.h"
0019 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0020 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0021 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0022 #include "FWCore/ServiceRegistry/interface/Service.h"
0023 #include "FWCore/Utilities/interface/InputTag.h"
0024 #include "FWCore/Utilities/interface/transform.h"
0025
0026 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0027 #include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h"
0028
0029 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0030 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0031
0032 #include "CLHEP/Units/GlobalSystemOfUnits.h"
0033 #include "CLHEP/Units/GlobalPhysicalConstants.h"
0034
0035 #include "TH1D.h"
0036 #include "TH2D.h"
0037
0038 class HGCalSimHitStudy : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0039 public:
0040 struct hitsinfo {
0041 hitsinfo() {
0042 phi = eta = energy = time = 0.0;
0043 layer = 0;
0044 }
0045 double phi, eta, energy, time;
0046 int layer;
0047 };
0048
0049 explicit HGCalSimHitStudy(const edm::ParameterSet&);
0050 ~HGCalSimHitStudy() override = default;
0051 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0052
0053 protected:
0054 void beginJob() override;
0055 void beginRun(edm::Run const&, edm::EventSetup const&) override;
0056 void analyze(edm::Event const&, edm::EventSetup const&) override;
0057 void endRun(edm::Run const&, edm::EventSetup const&) override {}
0058
0059 private:
0060 void analyzeHits(int, const std::string&, const std::vector<PCaloHit>&);
0061
0062
0063 const std::vector<std::string> nameDetectors_, caloHitSources_;
0064 const double rmin_, rmax_, zmin_, zmax_;
0065 const double etamin_, etamax_;
0066 const int nbinR_, nbinZ_, nbinEta_, nLayers_, verbosity_;
0067 const bool ifNose_, ifLayer_;
0068 const std::vector<edm::ESGetToken<HGCalDDDConstants, IdealGeometryRecord> > tok_hgcGeom_;
0069 const std::vector<edm::EDGetTokenT<edm::PCaloHitContainer> > tok_hits_;
0070 std::vector<const HGCalDDDConstants*> hgcons_;
0071 std::vector<int> layers_, layerFront_;
0072
0073
0074 std::vector<TH2D*> h_RZ_, h_EtaPhi_, h_EtFiZp_, h_EtFiZm_, h_XY_;
0075 std::vector<TH1D*> h_E_, h_T_, h_LayerZp_, h_LayerZm_;
0076 std::vector<TH1D*> h_W1_, h_W2_, h_C1_, h_C2_, h_Ly_;
0077 };
0078
0079 HGCalSimHitStudy::HGCalSimHitStudy(const edm::ParameterSet& iConfig)
0080 : nameDetectors_(iConfig.getParameter<std::vector<std::string> >("detectorNames")),
0081 caloHitSources_(iConfig.getParameter<std::vector<std::string> >("caloHitSources")),
0082 rmin_(iConfig.getUntrackedParameter<double>("rMin", 0.0)),
0083 rmax_(iConfig.getUntrackedParameter<double>("rMax", 3000.0)),
0084 zmin_(iConfig.getUntrackedParameter<double>("zMin", 3000.0)),
0085 zmax_(iConfig.getUntrackedParameter<double>("zMax", 6000.0)),
0086 etamin_(iConfig.getUntrackedParameter<double>("etaMin", 1.0)),
0087 etamax_(iConfig.getUntrackedParameter<double>("etaMax", 3.0)),
0088 nbinR_(iConfig.getUntrackedParameter<int>("nBinR", 300)),
0089 nbinZ_(iConfig.getUntrackedParameter<int>("nBinZ", 300)),
0090 nbinEta_(iConfig.getUntrackedParameter<int>("nBinEta", 200)),
0091 nLayers_(iConfig.getUntrackedParameter<int>("layers", 50)),
0092 verbosity_(iConfig.getUntrackedParameter<int>("verbosity", 0)),
0093 ifNose_(iConfig.getUntrackedParameter<bool>("ifNose", false)),
0094 ifLayer_(iConfig.getUntrackedParameter<bool>("ifLayer", false)),
0095 tok_hgcGeom_{
0096 edm::vector_transform(nameDetectors_,
0097 [this](const std::string& name) {
0098 return esConsumes<HGCalDDDConstants, IdealGeometryRecord, edm::Transition::BeginRun>(
0099 edm::ESInputTag{"", name});
0100 })},
0101 tok_hits_{edm::vector_transform(caloHitSources_, [this](const std::string& source) {
0102 return consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits", source));
0103 })} {
0104 usesResource(TFileService::kSharedResource);
0105 }
0106
0107 void HGCalSimHitStudy::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0108 edm::ParameterSetDescription desc;
0109 std::vector<std::string> names = {"HGCalEESensitive", "HGCalHESiliconSensitive", "HGCalHEScintillatorSensitive"};
0110 std::vector<std::string> sources = {"HGCHitsEE", "HGCHitsHEfront", "HGCHitsHEback"};
0111 desc.add<std::vector<std::string> >("detectorNames", names);
0112 desc.add<std::vector<std::string> >("caloHitSources", sources);
0113 desc.addUntracked<double>("rMin", 0.0);
0114 desc.addUntracked<double>("rMax", 3000.0);
0115 desc.addUntracked<double>("zMin", 3000.0);
0116 desc.addUntracked<double>("zMax", 6000.0);
0117 desc.addUntracked<double>("etaMin", 1.0);
0118 desc.addUntracked<double>("etaMax", 3.0);
0119 desc.addUntracked<int>("nBinR", 300);
0120 desc.addUntracked<int>("nBinZ", 300);
0121 desc.addUntracked<int>("nBinEta", 200);
0122 desc.addUntracked<int>("layers", 50);
0123 desc.addUntracked<int>("verbosity", 0);
0124 desc.addUntracked<bool>("ifNose", false);
0125 desc.addUntracked<bool>("ifLayer", false);
0126 descriptions.add("hgcalSimHitStudy", desc);
0127 }
0128
0129 void HGCalSimHitStudy::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0130
0131 for (unsigned int k = 0; k < tok_hits_.size(); ++k) {
0132 const edm::Handle<edm::PCaloHitContainer>& theCaloHitContainers = iEvent.getHandle(tok_hits_[k]);
0133 if (theCaloHitContainers.isValid()) {
0134 if (verbosity_ > 0)
0135 edm::LogVerbatim("HGCalValidation") << " PcalohitItr = " << theCaloHitContainers->size();
0136 std::vector<PCaloHit> caloHits;
0137 caloHits.insert(caloHits.end(), theCaloHitContainers->begin(), theCaloHitContainers->end());
0138 analyzeHits(k, nameDetectors_[k], caloHits);
0139 } else if (verbosity_ > 0) {
0140 edm::LogVerbatim("HGCalValidation") << "PCaloHitContainer does not "
0141 << "exist for " << nameDetectors_[k];
0142 }
0143 }
0144 }
0145
0146 void HGCalSimHitStudy::analyzeHits(int ih, std::string const& name, std::vector<PCaloHit> const& hits) {
0147 if (verbosity_ > 0)
0148 edm::LogVerbatim("HGCalValidation") << name << " with " << hits.size() << " PcaloHit elements";
0149
0150 std::map<uint32_t, hitsinfo> map_hits;
0151 map_hits.clear();
0152
0153 unsigned int nused(0);
0154 for (auto const& hit : hits) {
0155 double energy = hit.energy();
0156 double time = hit.time();
0157 uint32_t id = hit.id();
0158 int cell, sector, sector2(0), layer, zside;
0159 int subdet(0), cell2(0), type(0);
0160 HepGeom::Point3D<float> gcoord;
0161 std::pair<float, float> xy;
0162 if (ifNose_) {
0163 HFNoseDetId detId = HFNoseDetId(id);
0164 subdet = detId.subdetId();
0165 cell = detId.cellU();
0166 cell2 = detId.cellV();
0167 sector = detId.waferU();
0168 sector2 = detId.waferV();
0169 type = detId.type();
0170 layer = detId.layer();
0171 zside = detId.zside();
0172 xy = hgcons_[ih]->locateCell(zside, layer, sector, sector2, cell, cell2, false, true, false, false);
0173 h_W2_[ih]->Fill(sector2);
0174 h_C2_[ih]->Fill(cell2);
0175 } else if (hgcons_[ih]->waferHexagon8()) {
0176 HGCSiliconDetId detId = HGCSiliconDetId(id);
0177 subdet = static_cast<int>(detId.det());
0178 cell = detId.cellU();
0179 cell2 = detId.cellV();
0180 sector = detId.waferU();
0181 sector2 = detId.waferV();
0182 type = detId.type();
0183 layer = detId.layer();
0184 zside = detId.zside();
0185 xy = hgcons_[ih]->locateCell(zside, layer, sector, sector2, cell, cell2, false, true, false, false);
0186 h_W2_[ih]->Fill(sector2);
0187 h_C2_[ih]->Fill(cell2);
0188 } else if (hgcons_[ih]->tileTrapezoid()) {
0189 HGCScintillatorDetId detId = HGCScintillatorDetId(id);
0190 subdet = static_cast<int>(detId.det());
0191 sector = detId.ieta();
0192 cell = detId.iphi();
0193 type = detId.type();
0194 layer = detId.layer();
0195 zside = detId.zside();
0196 xy = hgcons_[ih]->locateCellTrap(zside, layer, sector, cell, false, false);
0197 } else {
0198 edm::LogError("HGCalValidation") << "HGCalSimHitStudy: Wrong geometry mode " << hgcons_[ih]->geomMode();
0199 continue;
0200 }
0201 double zp = hgcons_[ih]->waferZ(layer, false);
0202 if (zside < 0)
0203 zp = -zp;
0204 double xp = (zp < 0) ? -xy.first : xy.first;
0205 gcoord = HepGeom::Point3D<float>(xp, xy.second, zp);
0206 if (verbosity_ > 2)
0207 edm::LogVerbatim("HGCalValidation")
0208 << "i/p " << subdet << ":" << zside << ":" << layer << ":" << sector << ":" << sector2 << ":" << cell << ":"
0209 << cell2 << " o/p " << xy.first << ":" << xy.second << ":" << zp;
0210 nused++;
0211 double tof = (gcoord.mag() * CLHEP::mm) / CLHEP::c_light;
0212 if (verbosity_ > 1)
0213 edm::LogVerbatim("HGCalValidation")
0214 << "Detector " << name << " zside = " << zside << " layer = " << layer << " type = " << type
0215 << " wafer = " << sector << ":" << sector2 << " cell = " << cell << ":" << cell2 << " positon = " << gcoord
0216 << " energy = " << energy << " time = " << time << ":" << tof;
0217 time -= tof;
0218 if (time < 0)
0219 time = 0;
0220 hitsinfo hinfo;
0221 if (map_hits.count(id) != 0) {
0222 hinfo = map_hits[id];
0223 } else {
0224 hinfo.layer = layer + layerFront_[ih];
0225 hinfo.phi = gcoord.getPhi();
0226 hinfo.eta = gcoord.getEta();
0227 hinfo.time = time;
0228 }
0229 hinfo.energy += energy;
0230 map_hits[id] = hinfo;
0231
0232
0233 h_RZ_[0]->Fill(std::abs(gcoord.z()), gcoord.rho());
0234 h_RZ_[ih + 1]->Fill(std::abs(gcoord.z()), gcoord.rho());
0235 if (ifLayer_) {
0236 if (hinfo.layer <= static_cast<int>(h_XY_.size()))
0237 h_XY_[hinfo.layer - 1]->Fill(gcoord.x(), gcoord.y());
0238 } else {
0239 h_EtaPhi_[0]->Fill(std::abs(hinfo.eta), hinfo.phi);
0240 h_EtaPhi_[ih + 1]->Fill(std::abs(hinfo.eta), hinfo.phi);
0241 }
0242 h_Ly_[ih]->Fill(layer);
0243 h_W1_[ih]->Fill(sector);
0244 h_C1_[ih]->Fill(cell);
0245 }
0246 if (verbosity_ > 0)
0247 edm::LogVerbatim("HGCalValidation") << name << " with " << map_hits.size() << ":" << nused << " detector elements"
0248 << " being hit";
0249
0250 for (auto const& hit : map_hits) {
0251 hitsinfo hinfo = hit.second;
0252 if (verbosity_ > 1)
0253 edm::LogVerbatim("HGCalValidation")
0254 << " ---------------------- eta = " << hinfo.eta << " phi = " << hinfo.phi << " layer = " << hinfo.layer
0255 << " E = " << hinfo.energy << " T = " << hinfo.time;
0256 h_E_[0]->Fill(hinfo.energy);
0257 h_E_[ih + 1]->Fill(hinfo.energy);
0258 h_T_[0]->Fill(hinfo.time);
0259 h_T_[ih + 1]->Fill(hinfo.time);
0260 if (hinfo.eta > 0) {
0261 if (!ifLayer_) {
0262 h_EtFiZp_[0]->Fill(std::abs(hinfo.eta), hinfo.phi, hinfo.energy);
0263 h_EtFiZp_[ih + 1]->Fill(std::abs(hinfo.eta), hinfo.phi, hinfo.energy);
0264 }
0265 h_LayerZp_[0]->Fill(hinfo.layer, hinfo.energy);
0266 h_LayerZp_[ih + 1]->Fill(hinfo.layer, hinfo.energy);
0267 } else {
0268 if (!ifLayer_) {
0269 h_EtFiZm_[0]->Fill(std::abs(hinfo.eta), hinfo.phi, hinfo.energy);
0270 h_EtFiZm_[ih + 1]->Fill(std::abs(hinfo.eta), hinfo.phi, hinfo.energy);
0271 }
0272 h_LayerZm_[0]->Fill(hinfo.layer, hinfo.energy);
0273 h_LayerZm_[ih + 1]->Fill(hinfo.layer, hinfo.energy);
0274 }
0275 }
0276 }
0277
0278
0279 void HGCalSimHitStudy::beginRun(const edm::Run&, const edm::EventSetup& iSetup) {
0280 for (unsigned int k = 0; k < nameDetectors_.size(); ++k) {
0281 hgcons_.emplace_back(&iSetup.getData(tok_hgcGeom_[k]));
0282 layers_.emplace_back(hgcons_.back()->layers(false));
0283 layerFront_.emplace_back(hgcons_.back()->firstLayer());
0284 if (verbosity_ > 0)
0285 edm::LogVerbatim("HGCalValidation") << nameDetectors_[k] << " defined with " << layers_.back() << " Layers with "
0286 << layerFront_.back() << " in front";
0287 }
0288 }
0289
0290 void HGCalSimHitStudy::beginJob() {
0291 edm::Service<TFileService> fs;
0292
0293 std::ostringstream name, title;
0294 for (unsigned int ih = 0; ih <= nameDetectors_.size(); ++ih) {
0295 name.str("");
0296 title.str("");
0297 if (ih == 0) {
0298 name << "RZ_AllDetectors";
0299 title << "R vs Z for All Detectors";
0300 } else {
0301 name << "RZ_" << nameDetectors_[ih - 1];
0302 title << "R vs Z for " << nameDetectors_[ih - 1];
0303 }
0304 h_RZ_.emplace_back(
0305 fs->make<TH2D>(name.str().c_str(), title.str().c_str(), nbinZ_, zmin_, zmax_, nbinR_, rmin_, rmax_));
0306 if (ifLayer_) {
0307 if (ih == 0) {
0308 for (int ly = 0; ly < nLayers_; ++ly) {
0309 name.str("");
0310 title.str("");
0311 name << "XY_L" << (ly + 1);
0312 title << "Y vs X at Layer " << (ly + 1);
0313 h_XY_.emplace_back(
0314 fs->make<TH2D>(name.str().c_str(), title.str().c_str(), nbinR_, -rmax_, rmax_, nbinR_, -rmax_, rmax_));
0315 }
0316 }
0317 } else {
0318 name.str("");
0319 title.str("");
0320 if (ih == 0) {
0321 name << "EtaPhi_AllDetectors";
0322 title << "#phi vs #eta for All Detectors";
0323 } else {
0324 name << "EtaPhi_" << nameDetectors_[ih - 1];
0325 title << "#phi vs #eta for " << nameDetectors_[ih - 1];
0326 }
0327 h_EtaPhi_.emplace_back(
0328 fs->make<TH2D>(name.str().c_str(), title.str().c_str(), nbinEta_, etamin_, etamax_, 200, -M_PI, M_PI));
0329 name.str("");
0330 title.str("");
0331 if (ih == 0) {
0332 name << "EtFiZp_AllDetectors";
0333 title << "#phi vs #eta (+z) for All Detectors";
0334 } else {
0335 name << "EtFiZp_" << nameDetectors_[ih - 1];
0336 title << "#phi vs #eta (+z) for " << nameDetectors_[ih - 1];
0337 }
0338 h_EtFiZp_.emplace_back(
0339 fs->make<TH2D>(name.str().c_str(), title.str().c_str(), nbinEta_, etamin_, etamax_, 200, -M_PI, M_PI));
0340 name.str("");
0341 title.str("");
0342 if (ih == 0) {
0343 name << "EtFiZm_AllDetectors";
0344 title << "#phi vs #eta (-z) for All Detectors";
0345 } else {
0346 name << "EtFiZm_" << nameDetectors_[ih - 1];
0347 title << "#phi vs #eta (-z) for " << nameDetectors_[ih - 1];
0348 }
0349 h_EtFiZm_.emplace_back(
0350 fs->make<TH2D>(name.str().c_str(), title.str().c_str(), nbinEta_, etamin_, etamax_, 200, -M_PI, M_PI));
0351 }
0352 name.str("");
0353 title.str("");
0354 if (ih == 0) {
0355 name << "LayerZp_AllDetectors";
0356 title << "Energy vs Layer (+z) for All Detectors";
0357 } else {
0358 name << "LayerZp_" << nameDetectors_[ih - 1];
0359 title << "Energy vs Layer (+z) for " << nameDetectors_[ih - 1];
0360 }
0361 h_LayerZp_.emplace_back(fs->make<TH1D>(name.str().c_str(), title.str().c_str(), 60, 0.0, 60.0));
0362 name.str("");
0363 title.str("");
0364 if (ih == 0) {
0365 name << "LayerZm_AllDetectors";
0366 title << "Energy vs Layer (-z) for All Detectors";
0367 } else {
0368 name << "LayerZm_" << nameDetectors_[ih - 1];
0369 title << "Energy vs Layer (-z) for " << nameDetectors_[ih - 1];
0370 }
0371 h_LayerZm_.emplace_back(fs->make<TH1D>(name.str().c_str(), title.str().c_str(), 60, 0.0, 60.0));
0372
0373 name.str("");
0374 title.str("");
0375 if (ih == 0) {
0376 name << "E_AllDetectors";
0377 title << "Energy Deposit for All Detectors";
0378 } else {
0379 name << "E_" << nameDetectors_[ih - 1];
0380 title << "Energy Deposit for " << nameDetectors_[ih - 1];
0381 }
0382 h_E_.emplace_back(fs->make<TH1D>(name.str().c_str(), title.str().c_str(), 1000, 0.0, 1.0));
0383
0384 name.str("");
0385 title.str("");
0386 if (ih == 0) {
0387 name << "T_AllDetectors";
0388 title << "Time for All Detectors";
0389 } else {
0390 name << "T_" << nameDetectors_[ih - 1];
0391 title << "Time for " << nameDetectors_[ih - 1];
0392 }
0393 h_T_.emplace_back(fs->make<TH1D>(name.str().c_str(), title.str().c_str(), 1000, 0.0, 200.0));
0394 }
0395
0396 for (unsigned int ih = 0; ih < nameDetectors_.size(); ++ih) {
0397 name.str("");
0398 title.str("");
0399 name << "LY_" << nameDetectors_[ih];
0400 title << "Layer number for " << nameDetectors_[ih];
0401 h_Ly_.emplace_back(fs->make<TH1D>(name.str().c_str(), title.str().c_str(), 200, 0, 100));
0402 if (nameDetectors_[ih] == "HGCalHEScintillatorSensitive") {
0403 name.str("");
0404 title.str("");
0405 name << "IR_" << nameDetectors_[ih];
0406 title << "Radius index for " << nameDetectors_[ih];
0407 h_W1_.emplace_back(fs->make<TH1D>(name.str().c_str(), title.str().c_str(), 200, -50, 50));
0408 name.str("");
0409 title.str("");
0410 name << "FI_" << nameDetectors_[ih];
0411 title << "#phi index for " << nameDetectors_[ih];
0412 h_C1_.emplace_back(fs->make<TH1D>(name.str().c_str(), title.str().c_str(), 720, 0, 360));
0413 } else {
0414 name.str("");
0415 title.str("");
0416 name << "WU_" << nameDetectors_[ih];
0417 title << "u index of wafers for " << nameDetectors_[ih];
0418 h_W1_.emplace_back(fs->make<TH1D>(name.str().c_str(), title.str().c_str(), 200, -50, 50));
0419 name.str("");
0420 title.str("");
0421 name << "WV_" << nameDetectors_[ih];
0422 title << "v index of wafers for " << nameDetectors_[ih];
0423 h_W2_.emplace_back(fs->make<TH1D>(name.str().c_str(), title.str().c_str(), 100, -50, 50));
0424 name.str("");
0425 title.str("");
0426 name << "CU_" << nameDetectors_[ih];
0427 title << "u index of cells for " << nameDetectors_[ih];
0428 h_C1_.emplace_back(fs->make<TH1D>(name.str().c_str(), title.str().c_str(), 100, 0, 50));
0429 name.str("");
0430 title.str("");
0431 name << "CV_" << nameDetectors_[ih];
0432 title << "v index of cells for " << nameDetectors_[ih];
0433 h_C2_.emplace_back(fs->make<TH1D>(name.str().c_str(), title.str().c_str(), 100, 0, 50));
0434 }
0435 }
0436 }
0437
0438 #include "FWCore/Framework/interface/MakerMacros.h"
0439
0440 DEFINE_FWK_MODULE(HGCalSimHitStudy);
