Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:15:09

0001 // -*- C++ -*-
0002 //
0003 // Package:    HGCalValidityTester
0004 // Class:      HGCalValidityTester
0005 //
0006 /**\class HGCalValidityTester HGCalValidityTester.cc
0007  test/HGCalValidityTester.cc
0008 
0009  Description: <one line class summary>
0010 
0011  Implementation:
0012      <Notes on implementation>
0013 */
0014 //
0015 // Original Author:  Sunanda Banerjee
0016 //         Created:  Mon 2022/10/22
0017 //
0018 //
0019 
0020 // system include files
0021 #include <fstream>
0022 #include <iostream>
0023 #include <sstream>
0024 #include <string>
0025 #include <vector>
0026 
0027 // user include files
0028 #include "FWCore/Framework/interface/Frameworkfwd.h"
0029 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0030 #include "FWCore/Framework/interface/Event.h"
0031 #include "FWCore/Framework/interface/EventSetup.h"
0032 #include "FWCore/Framework/interface/MakerMacros.h"
0033 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0034 #include "FWCore/ParameterSet/interface/FileInPath.h"
0035 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0036 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0037 #include "FWCore/Utilities/interface/transform.h"
0038 
0039 #include "DataFormats/DetId/interface/DetId.h"
0040 #include "DataFormats/ForwardDetId/interface/HGCSiliconDetId.h"
0041 #include "DataFormats/ForwardDetId/interface/HGCScintillatorDetId.h"
0042 #include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h"
0043 #include "Geometry/HGCalCommonData/interface/HGCalGeomUtils.h"
0044 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0045 
0046 class HGCalValidityTester : public edm::one::EDAnalyzer<edm::one::WatchRuns> {
0047 public:
0048   explicit HGCalValidityTester(const edm::ParameterSet &);
0049   ~HGCalValidityTester() override = default;
0050   static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0051 
0052   void beginJob() override {}
0053   void beginRun(edm::Run const &, edm::EventSetup const &) override;
0054   void analyze(edm::Event const &iEvent, edm::EventSetup const &) override {}
0055   void endRun(edm::Run const &, edm::EventSetup const &) override {}
0056   void endJob() override {}
0057 
0058 private:
0059   const std::vector<std::string> nameDetectors_;
0060   const std::string fileName_;
0061   const std::vector<edm::ESGetToken<HGCalDDDConstants, IdealGeometryRecord>> tok_hgcal_;
0062   std::vector<const HGCalDDDConstants *> hgcCons_;
0063   std::vector<std::pair<DetId, uint32_t>> detIds_;
0064 };
0065 
0066 HGCalValidityTester::HGCalValidityTester(const edm::ParameterSet &iC)
0067     : nameDetectors_(iC.getParameter<std::vector<std::string>>("nameDetectors")),
0068       fileName_(iC.getParameter<std::string>("fileName")),
0069       tok_hgcal_{edm::vector_transform(nameDetectors_, [this](const std::string &name) {
0070         return esConsumes<HGCalDDDConstants, IdealGeometryRecord, edm::Transition::BeginRun>(edm::ESInputTag{"", name});
0071       })} {
0072   std::ostringstream st1;
0073   for (const auto &name : nameDetectors_)
0074     st1 << " : " << name;
0075   edm::LogVerbatim("HGCGeom") << "Test validity of cells for " << nameDetectors_.size() << " detectors" << st1.str()
0076                               << " with inputs from " << fileName_;
0077   if (!fileName_.empty()) {
0078     edm::FileInPath filetmp("Geometry/HGCalCommonData/data/" + fileName_);
0079     std::string fileName = filetmp.fullPath();
0080     std::ifstream fInput(fileName.c_str());
0081     if (!fInput.good()) {
0082       edm::LogVerbatim("HGCGeom") << "Cannot open file " << fileName;
0083     } else {
0084       char buffer[80];
0085       const std::vector<DetId::Detector> dets = {DetId::HGCalEE, DetId::HGCalHSi, DetId::HGCalHSc};
0086       while (fInput.getline(buffer, 80)) {
0087         std::vector<std::string> items = HGCalGeomUtils::splitString(std::string(buffer));
0088         if (items.size() == 8) {
0089           DetId::Detector det = static_cast<DetId::Detector>(std::atoi(items[0].c_str()));
0090           auto itr = std::find(dets.begin(), dets.end(), det);
0091           if (itr != dets.end()) {
0092             uint32_t pos = static_cast<uint32_t>(itr - dets.begin());
0093             DetId id(0);
0094             if ((det == DetId::HGCalEE) || (det == DetId::HGCalHSi)) {
0095               int type = std::atoi(items[1].c_str());
0096               int zside = std::atoi(items[2].c_str());
0097               int layer = std::atoi(items[3].c_str());
0098               int waferU = std::atoi(items[4].c_str());
0099               int waferV = std::atoi(items[5].c_str());
0100               int cellU = std::atoi(items[6].c_str());
0101               int cellV = std::atoi(items[7].c_str());
0102               id = static_cast<DetId>(HGCSiliconDetId(det, zside, type, layer, waferU, waferV, cellU, cellV));
0103             } else if (det == DetId::HGCalHSc) {
0104               bool trig = (std::atoi(items[1].c_str()) > 0);
0105               int type = std::atoi(items[2].c_str());
0106               int zside = std::atoi(items[3].c_str());
0107               int sipm = std::atoi(items[4].c_str());
0108               int layer = std::atoi(items[5].c_str());
0109               int ring = std::atoi(items[6].c_str());
0110               int iphi = std::atoi(items[7].c_str());
0111               id = static_cast<DetId>(HGCScintillatorDetId(type, layer, zside * ring, iphi, trig, sipm));
0112             }
0113             if (id.rawId() != 0)
0114               detIds_.emplace_back(id, pos);
0115           }
0116         }
0117       }
0118       fInput.close();
0119     }
0120   }
0121   edm::LogVerbatim("HGCGeom") << "Reads " << detIds_.size() << " ID's from " << fileName_;
0122   for (unsigned int k = 0; k < detIds_.size(); ++k) {
0123     if ((detIds_[k].first).det() == DetId::HGCalHSc)
0124       edm::LogVerbatim("HGCalGeom") << "[" << k << "] " << HGCScintillatorDetId(detIds_[k].first) << " from DDConstant "
0125                                     << (detIds_[k].second);
0126     else
0127       edm::LogVerbatim("HGCalGeom") << "[" << k << "] " << HGCSiliconDetId(detIds_[k].first) << " from DDConstant "
0128                                     << (detIds_[k].second);
0129   }
0130 }
0131 
0132 void HGCalValidityTester::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0133   std::vector<std::string> names = {"HGCalEESensitive", "HGCalHESiliconSensitive", "HGCalHEScintillatorSensitive"};
0134   edm::ParameterSetDescription desc;
0135   desc.add<std::vector<std::string>>("nameDetectors", names);
0136   desc.add<std::string>("fileName", "missD88.txt");
0137   descriptions.add("hgcalValidityTester", desc);
0138 }
0139 
0140 // ------------ method called to produce the data  ------------
0141 void HGCalValidityTester::beginRun(edm::Run const &iRun, edm::EventSetup const &iSetup) {
0142   //initiating hgc Geometry
0143   std::vector<std::string> names = {"HGCalEESensitive", "HGCalHESiliconSensitive", "HGCalHEScintillatorSensitive"};
0144   std::vector<DetId::Detector> dets = {DetId::HGCalEE, DetId::HGCalHSi, DetId::HGCalHSc};
0145   std::map<DetId::Detector, uint32_t> detMap;
0146   for (uint32_t i = 0; i < nameDetectors_.size(); i++) {
0147     edm::LogVerbatim("HGCGeom") << "Tries to initialize HGCalGeometry and HGCalDDDConstants for " << i << ":"
0148                                 << nameDetectors_[i];
0149     const edm::ESHandle<HGCalDDDConstants> &hgcCons = iSetup.getHandle(tok_hgcal_[i]);
0150     if (hgcCons.isValid()) {
0151       hgcCons_.push_back(hgcCons.product());
0152     } else {
0153       edm::LogWarning("HGCGeom") << "Cannot initiate HGCalDDDConstants for " << nameDetectors_[i] << std::endl;
0154     }
0155     auto ii = std::find(names.begin(), names.end(), nameDetectors_[i]);
0156     if (ii != names.end()) {
0157       uint32_t k = static_cast<uint32_t>(ii - names.begin());
0158       detMap[dets[k]] = i;
0159     }
0160   }
0161   edm::LogVerbatim("HGCGeom") << "Loaded HGCalDDConstants for " << detMap.size() << " detectors";
0162 
0163   for (auto itr = detMap.begin(); itr != detMap.end(); ++itr)
0164     edm::LogVerbatim("HGCGeom") << "[" << itr->second << "]: " << nameDetectors_[itr->second] << " for Detector "
0165                                 << itr->first;
0166 
0167   for (unsigned int k = 0; k < detIds_.size(); ++k) {
0168     std::ostringstream st1;
0169     std::pair<float, float> xy;
0170     bool valid(false);
0171     int layer(0), zside(1);
0172     const HGCalDDDConstants *cons = hgcCons_[detMap[(detIds_[k].first).det()]];
0173     if ((detIds_[k].first).det() == DetId::HGCalHSc) {
0174       HGCScintillatorDetId id(detIds_[k].first);
0175       layer = id.layer();
0176       zside = id.zside();
0177       xy = cons->locateCellTrap(zside, layer, id.ring(), id.iphi(), true, false);
0178       double z = zside * (cons->waferZ(layer, true));
0179       valid = cons->isValidTrap(zside, layer, id.ring(), id.iphi());
0180       auto cell = cons->assignCellTrap(zside * xy.first, xy.second, z, layer, true);
0181       HGCScintillatorDetId newId(cell[2], layer, zside * cell[0], cell[1], id.trigger(), id.sipm());
0182       st1 << "Old: " << id << " New: " << newId << " OK " << (id.rawId() == newId.rawId());
0183     } else {
0184       HGCSiliconDetId id(detIds_[k].first);
0185       layer = id.layer();
0186       zside = id.zside();
0187       xy = cons->locateCell(zside, layer, id.waferU(), id.waferV(), id.cellU(), id.cellV(), true, true, false, false);
0188       valid = cons->isValidHex8(layer, id.waferU(), id.waferV(), id.cellU(), id.cellV(), false);
0189       auto cell = cons->assignCellHex(zside * xy.first, xy.second, zside, layer, true, false, false);
0190       HGCSiliconDetId newId(id.det(), id.zside(), cell[2], id.layer(), cell[0], cell[1], cell[3], cell[4]);
0191       st1 << "Old: " << id << " New: " << newId << " OK " << (id.rawId() == newId.rawId());
0192     }
0193     double r = std::sqrt(xy.first * xy.first + xy.second * xy.second);
0194     double z = zside * (cons->waferZ(layer, true));
0195     auto range = cons->getRangeR(layer, true);
0196     edm::LogVerbatim("HGCalMiss") << "Hit[" << k << "] " << st1.str() << " Position (" << xy.first << ", " << xy.second
0197                                   << ", " << z << ") Valid " << valid << " R " << r << " (" << range.first << ":"
0198                                   << range.second << ")";
0199   }
0200 }
0201 
0202 // define this as a plug-in
0203 DEFINE_FWK_MODULE(HGCalValidityTester);