Calo/plugins/HGCalTestGuardRing.cc

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/MessageLogger/interface/MessageLogger.h"
0009 #include "FWCore/ParameterSet/interface/FileInPath.h"
0010 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0011 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0012 #include "FWCore/Utilities/interface/InputTag.h"
0013 #include "FWCore/Utilities/interface/Exception.h"
0014
0015 #include "DataFormats/ForwardDetId/interface/HGCSiliconDetId.h"
0016 #include "DataFormats/ForwardDetId/interface/HGCScintillatorDetId.h"
0017
0018 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0019 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0020 #include "SimG4CMS/Calo/interface/CaloSimUtils.h"
0021 #include "SimG4CMS/Calo/interface/HGCGuardRing.h"
0022
0023 #include "Geometry/CaloTopology/interface/HGCalTopology.h"
0024 #include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
0025 #include "Geometry/HGCalCommonData/interface/HGCalCell.h"
0026 #include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h"
0027 #include "Geometry/HGCalCommonData/interface/HGCalParameters.h"
0028 #include "Geometry/HGCalCommonData/interface/HGCalTypes.h"
0029 #include "Geometry/HGCalCommonData/interface/HGCalWaferMask.h"
0030 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0031
0032 #include <fstream>
0033 #include <sstream>
0034 #include <string>
0035 #include <map>
0036
0037 class HGCalTestGuardRing : public edm::one::EDAnalyzer<> {
0038 public:
0039   HGCalTestGuardRing(const edm::ParameterSet& ps);
0040   ~HGCalTestGuardRing() override = default;
0041   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0042
0043 protected:
0044   void analyze(edm::Event const&, edm::EventSetup const&) override;
0045   void beginJob() override {}
0046   void endJob() override {}
0047
0048 private:
0049   const std::string nameSense_, waferFile_;
0050   const double guardRingOffset_;
0051   const edm::ESGetToken<HGCalGeometry, IdealGeometryRecord> geomToken_;
0052   std::map<HGCSiliconDetId, int> waferID_;
0053 };
0054
0055 HGCalTestGuardRing::HGCalTestGuardRing(const edm::ParameterSet& ps)
0056     : nameSense_(ps.getParameter<std::string>("nameSense")),
0057       waferFile_(ps.getParameter<std::string>("waferFile")),
0058       guardRingOffset_(ps.getParameter<double>("guardRingOffset")),
0059       geomToken_(esConsumes<HGCalGeometry, IdealGeometryRecord>(edm::ESInputTag{"", nameSense_})) {
0060   DetId::Detector det = (nameSense_ != "HGCalHESiliconSensitive") ? DetId::HGCalEE : DetId::HGCalHSi;
0061   edm::LogVerbatim("HGCalSim") << "Test Guard Ring Offset " << guardRingOffset_ << " for " << nameSense_ << ":" << det
0062                                << " for wafers read from file " << waferFile_;
0063   if (!waferFile_.empty()) {
0064     std::string thick[4] = {"h120", "l200", "l300", "h200"};
0065     int addType[4] = {HGCalTypes::WaferHD120, HGCalTypes::WaferLD200, HGCalTypes::WaferLD300, HGCalTypes::WaferHD200};
0066     const int partTypeH[6] = {HGCalTypes::WaferFull,
0067                               HGCalTypes::WaferHalf2,
0068                               HGCalTypes::WaferChopTwoM,
0069                               HGCalTypes::WaferSemi2,
0070                               HGCalTypes::WaferSemi2,
0071                               HGCalTypes::WaferFive2};
0072     const int partTypeL[7] = {HGCalTypes::WaferFull,
0073                               HGCalTypes::WaferHalf,
0074                               HGCalTypes::WaferHalf,
0075                               HGCalTypes::WaferSemi,
0076                               HGCalTypes::WaferSemi,
0077                               HGCalTypes::WaferFive,
0078                               HGCalTypes::WaferThree};
0079     edm::FileInPath filetmp("SimG4CMS/Calo/data/" + waferFile_);
0080     const std::string& fileName = filetmp.fullPath();
0081     std::ifstream fInput(fileName.c_str());
0082     if (!fInput.good()) {
0083       edm::LogVerbatim("HGCalSim") << "Cannot open file " << fileName;
0084     } else {
0085       char buffer[80];
0086       while (fInput.getline(buffer, 80)) {
0087         std::vector<std::string> items = CaloSimUtils::splitString(std::string(buffer));
0088         if (items.size() > 6) {
0089           int layer = std::atoi(items[0].c_str());
0090           int waferU = std::atoi(items[4].c_str());
0091           int waferV = std::atoi(items[5].c_str());
0092           int thck = static_cast<int>(std::find(thick, thick + 4, items[2]) - thick);
0093           int type = (thck < 4) ? addType[thck] : 0;
0094           HGCSiliconDetId id(det, -1, type, layer, waferU, waferV, 0, 0);
0095           int orient = std::atoi(items[5].c_str());
0096           int part = std::atoi(items[1].c_str());
0097           if (part >= 0) {
0098             if ((type == HGCalTypes::WaferHD120) || (type == HGCalTypes::WaferHD200))
0099               part = partTypeH[part];
0100             else
0101               part = partTypeL[part];
0102           }
0103           waferID_[id] = orient * 100 + part;
0104 #ifdef EDM_ML_DEBUG
0105           edm::LogVerbatim("HGCalSim") << "HGCalTestGuardRing::Reads " << id << " Orientation:Partial " << orient << ":"
0106                                        << part;
0107 #endif
0108         }
0109       }
0110       edm::LogVerbatim("HGCalSim") << "HGCalTestGuardRing::Reads in " << waferID_.size() << " wafer information from "
0111                                    << fileName;
0112       fInput.close();
0113     }
0114   }
0115 }
0116
0117 void HGCalTestGuardRing::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0118   edm::ParameterSetDescription desc;
0119   desc.add<std::string>("nameSense", "HGCalEESensitive");
0120   desc.add<std::string>("waferFile", "testWafersEE.txt");
0121   desc.add<double>("guardRingOffset", 1.0);
0122   descriptions.add("hgcalTestGuardRingEE", desc);
0123 }
0124
0125 void HGCalTestGuardRing::analyze(const edm::Event& e, const edm::EventSetup& iS) {
0126   // get HGCalGeometry
0127   const HGCalGeometry* geom = &iS.getData(geomToken_);
0128   const HGCalDDDConstants& hgc = geom->topology().dddConstants();
0129   double waferSize = hgc.waferSize(false);
0130   HGCalCell wafer(waferSize, hgc.getUVMax(0), hgc.getUVMax(1));
0131   const bool v17OrLess = hgc.v17OrLess();
0132
0133   // get the hit collection
0134   edm::LogVerbatim("HGCalSim") << "HGCalTestGuardRing: Wafer Szie " << waferSize << " v17OrLess " << v17OrLess;
0135
0136   // Loop over all IDs
0137   int all(0), allSi(0), good(0);
0138   for (std::map<HGCSiliconDetId, int>::const_iterator itr = waferID_.begin(); itr != waferID_.end(); ++itr) {
0139     HGCSiliconDetId id = itr->first;
0140     ++all;
0141     if ((id.det() == DetId::HGCalEE) || (id.det() == DetId::HGCalHSi)) {
0142       ++allSi;
0143       if (((id.det() == DetId::HGCalEE) && (nameSense_ == "HGCalEESensitive")) ||
0144           ((id.det() == DetId::HGCalHSi) && (nameSense_ == "HGCalHESiliconSensitive"))) {
0145         int partial = ((itr->second) % 100);
0146         int orient = (((itr->second) / 100) % 100);
0147         int type = id.type();
0148         int nCells = hgc.getUVMax(type);
0149         for (int u = 0; u < 2 * nCells; ++u) {
0150           for (int v = 0; v < 2 * nCells; ++v) {
0151             if (((v - u) < nCells) && ((u - v) <= nCells)) {
0152               HGCSiliconDetId hid(id.det(), id.zside(), id.type(), id.layer(), id.waferU(), id.waferV(), u, v);
0153               bool valid = (geom->topology()).valid(static_cast<DetId>(hid));
0154               if (valid) {
0155                 ++good;
0156                 int placeIndex = wafer.cellPlacementIndex(1, HGCalTypes::waferFrontBack(0), orient);
0157                 std::pair<double, double> xy = wafer.cellUV2XY1(u, v, placeIndex, type);
0158                 std::vector<std::pair<double, double> > wxy1 =
0159                     HGCalWaferMask::waferXY(partial, orient, -1, waferSize, 0.0, 0.0, 0.0, v17OrLess);
0160                 bool check1 = HGCGuardRing::insidePolygon(xy.first, xy.second, wxy1);
0161                 std::ostringstream st1;
0162                 for (unsigned int k1 = 0; k1 < wxy1.size(); ++k1)
0163                   st1 << " (" << wxy1[k1].first << ", " << wxy1[k1].second << ")";
0164                 edm::LogVerbatim("HGCSim")
0165                     << "First " << hid << " Type:Partial:Orient:Place " << type << ":" << partial << ":" << orient
0166                     << ":" << placeIndex << " Boundary with " << wxy1.size() << " points: " << st1.str() << " check "
0167                     << check1 << " for (" << xy.first << ", " << xy.second << ")";
0168
0169                 std::vector<std::pair<double, double> > wxy2 =
0170                     HGCalWaferMask::waferXY(partial, orient, -1, waferSize, guardRingOffset_, 0.0, 0.0, v17OrLess);
0171                 bool check2 = HGCGuardRing::insidePolygon(xy.first, xy.second, wxy2);
0172                 std::ostringstream st2;
0173                 for (unsigned int k1 = 0; k1 < wxy2.size(); ++k1)
0174                   st2 << " (" << wxy2[k1].first << ", " << wxy2[k1].second << ")";
0175                 edm::LogVerbatim("HGCSim") << "Second Offset " << guardRingOffset_ << " Boundary with " << wxy2.size()
0176                                            << " points: " << st2.str() << " check " << check2 << " for (" << xy.first
0177                                            << ", " << xy.second << ")";
0178               }
0179             }
0180           }
0181         }
0182       }
0183     }
0184   }
0185   edm::LogVerbatim("HGCalSim") << "Total hits = " << all << " Good Silicon DetIds = " << allSi << ":" << good;
0186 }
0187
0188 //define this as a plug-in
0189 DEFINE_FWK_MODULE(HGCalTestGuardRing);