Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-12-20 03:14:12 

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 
0022 #include "Geometry/CaloTopology/interface/HGCalTopology.h"
0023 #include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
0024 #include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h"
0025 #include "Geometry/HGCalCommonData/interface/HGCalParameters.h"
0026 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0027 
0028 #include <fstream>
0029 #include <sstream>
0030 #include <string>
0031 #include <vector>
0032 
0033 class HGCalTestPartialWaferHits : public edm::one::EDAnalyzer<> {
0034 public:
0035   HGCalTestPartialWaferHits(const edm::ParameterSet& ps);
0036   ~HGCalTestPartialWaferHits() override = default;
0037   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0038 
0039 protected:
0040   void analyze(edm::Event const&, edm::EventSetup const&) override;
0041   void beginJob() override {}
0042   void endJob() override {}
0043 
0044 private:
0045   const std::string g4Label_, caloHitSource_, nameSense_, missingFile_;
0046   const edm::EDGetTokenT<edm::PCaloHitContainer> tok_calo_;
0047   const edm::ESGetToken<HGCalGeometry, IdealGeometryRecord> geomToken_;
0048   std::vector<int> wafers_;
0049   std::vector<int> dumpDets_;
0050 };
0051 
0052 HGCalTestPartialWaferHits::HGCalTestPartialWaferHits(const edm::ParameterSet& ps)
0053     : g4Label_(ps.getParameter<std::string>("moduleLabel")),
0054       caloHitSource_(ps.getParameter<std::string>("caloHitSource")),
0055       nameSense_(ps.getParameter<std::string>("nameSense")),
0056       missingFile_(ps.getParameter<std::string>("missingFile")),
0057       tok_calo_(consumes<edm::PCaloHitContainer>(edm::InputTag(g4Label_, caloHitSource_))),
0058       geomToken_(esConsumes<HGCalGeometry, IdealGeometryRecord>(edm::ESInputTag{"", nameSense_})) {
0059   edm::LogVerbatim("HGCalSim") << "Test Hit ID using SimHits for " << nameSense_ << " with module Label: " << g4Label_
0060                                << "   Hits: " << caloHitSource_ << " Missing Wafer file " << missingFile_;
0061   if (!missingFile_.empty()) {
0062     edm::FileInPath filetmp("SimG4CMS/Calo/data/" + missingFile_);
0063     const std::string& fileName = filetmp.fullPath();
0064     std::ifstream fInput(fileName.c_str());
0065     if (!fInput.good()) {
0066       edm::LogVerbatim("HGCalSim") << "Cannot open file " << fileName;
0067     } else {
0068       char buffer[80];
0069       while (fInput.getline(buffer, 80)) {
0070         std::vector<std::string> items = CaloSimUtils::splitString(std::string(buffer));
0071         if (items.size() > 2) {
0072           int layer = std::atoi(items[0].c_str());
0073           int waferU = std::atoi(items[1].c_str());
0074           int waferV = std::atoi(items[2].c_str());
0075           wafers_.emplace_back(HGCalWaferIndex::waferIndex(layer, waferU, waferV, false));
0076         } else if (items.size() == 1) {
0077           int dumpdet = std::atoi(items[0].c_str());
0078           dumpDets_.emplace_back(dumpdet);
0079           edm::LogVerbatim("HGCalSim") << nameSense_ << " Dump detector " << dumpdet;
0080         }
0081       }
0082       edm::LogVerbatim("HGCalSim") << "HGCalTestPartialWaferHits::Reads in " << wafers_.size() << ":"
0083                                    << dumpDets_.size() << " wafer|detector information from " << fileName;
0084       fInput.close();
0085     }
0086   }
0087 }
0088 
0089 void HGCalTestPartialWaferHits::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0090   edm::ParameterSetDescription desc;
0091   desc.add<std::string>("moduleLabel", "g4SimHits");
0092   desc.add<std::string>("caloHitSource", "HGCHitsEE");
0093   desc.add<std::string>("nameSense", "HGCalEESensitive");
0094   desc.add<std::string>("missingFile", "");
0095   descriptions.add("hgcalHitPartialEE", desc);
0096 }
0097 
0098 void HGCalTestPartialWaferHits::analyze(const edm::Event& e, const edm::EventSetup& iS) {
0099   // get HGCalGeometry
0100   const HGCalGeometry* geom = &iS.getData(geomToken_);
0101   const HGCalDDDConstants& hgc = geom->topology().dddConstants();
0102   int firstLayer = hgc.getLayerOffset();
0103   // get the hit collection
0104   const edm::Handle<edm::PCaloHitContainer>& hitsCalo = e.getHandle(tok_calo_);
0105   bool getHits = (hitsCalo.isValid());
0106   uint32_t nhits = (getHits) ? hitsCalo->size() : 0;
0107   uint32_t good(0), allSi(0), all(0), allSc(0), bad(0);
0108   constexpr double tol = 2.0;
0109   edm::LogVerbatim("HGCalSim") << "HGCalTestPartialWaferHits: Input flags Hits " << getHits << " with " << nhits
0110                                << " hits: Layer Offset " << firstLayer;
0111 
0112   if (getHits) {
0113     std::vector<PCaloHit> hits;
0114     hits.insert(hits.end(), hitsCalo->begin(), hitsCalo->end());
0115     if (!hits.empty()) {
0116       // Loop over all hits
0117       for (auto hit : hits) {
0118         ++all;
0119         DetId id(hit.id());
0120         bool valid = (geom->topology()).valid(id);
0121         std::ostringstream st1;
0122         if ((id.det() == DetId::HGCalEE) || (id.det() == DetId::HGCalHSi)) {
0123           ++allSi;
0124           HGCSiliconDetId hid(id);
0125           st1 << hid;
0126           if (((id.det() == DetId::HGCalEE) && (nameSense_ == "HGCalEESensitive")) ||
0127               ((id.det() == DetId::HGCalHSi) && (nameSense_ == "HGCalHESiliconSensitive"))) {
0128             std::string_view pid =
0129                 ((hgc.cassetteShiftSilicon(hid.zside(), hid.layer(), hid.waferU(), hid.waferV())) ? "HGCSim"
0130                                                                                                   : "HGCalSim");
0131             const auto& info = hgc.waferInfo(hid.layer(), hid.waferU(), hid.waferV());
0132             if (!valid)
0133               st1 << " Wafer Type:Part:Orient:Cassette " << info.type << ":" << info.part << ":" << info.orient << ":"
0134                   << info.cassette;
0135             bool toCheck(false);
0136             if (!wafers_.empty()) {
0137               int indx = HGCalWaferIndex::waferIndex(firstLayer + hid.layer(), hid.waferU(), hid.waferV(), false);
0138               if (std::find(wafers_.begin(), wafers_.end(), indx) != wafers_.end())
0139                 toCheck = true;
0140             } else if (!dumpDets_.empty()) {
0141               if ((std::find(dumpDets_.begin(), dumpDets_.end(), static_cast<int>(id.det())) != dumpDets_.end()) &&
0142                   (info.part != HGCalTypes::WaferFull))
0143                 toCheck = true;
0144             } else {
0145               // Only partial wafers
0146               toCheck = (info.part != HGCalTypes::WaferFull);
0147             }
0148             if (toCheck) {
0149               ++good;
0150               GlobalPoint pos = geom->getPosition(id);
0151               bool valid1 = geom->topology().valid(id);
0152               bool valid2 = hgc.isValidHex8(hid.layer(), hid.waferU(), hid.waferV(), hid.cellU(), hid.cellV(), false);
0153               auto xy = hgc.locateCell(hid, false);
0154               double xx = (hid.zside() > 0) ? xy.first : -xy.first;
0155               double dx = xx - pos.x();
0156               double dy = xy.second - pos.y();
0157               double diff = std::sqrt(dx * dx + dy * dy);
0158               if ((diff > tol) || (!valid1) || (!valid2))
0159                 pid = "HGCalError";
0160               edm::LogVerbatim(pid) << "Hit[" << all << ":" << allSi << ":" << good << "]" << hid
0161                                     << " Wafer Type:Part:Orient:Cassette " << info.type << ":" << info.part << ":"
0162                                     << info.orient << ":" << info.cassette << " at (" << pos.x() << ":" << xx << ":"
0163                                     << dx << ", " << pos.y() << ":" << xy.second << ":" << dy << ", " << pos.z()
0164                                     << ") Valid " << valid1 << ":" << valid2 << " Distance " << diff;
0165             }
0166           }
0167         } else if (id.det() == DetId::HGCalHSc) {
0168           ++allSc;
0169           HGCScintillatorDetId hid(id);
0170           st1 << hid;
0171           if ((id.det() == DetId::HGCalHSc) && (nameSense_ == "HGCalHEScintillatorSensitive")) {
0172             std::string_view pid =
0173                 ((hgc.cassetteShiftScintillator(hid.zside(), hid.layer(), hid.iphi())) ? "HGCSim" : "HGCalSim");
0174             GlobalPoint pos = geom->getPosition(id);
0175             bool valid1 = geom->topology().valid(id);
0176             bool valid2 = hgc.isValidTrap(hid.zside(), hid.layer(), hid.ring(), hid.iphi());
0177             if ((!valid1) || (!valid2))
0178               pid = "HGCalError";
0179             int cassette = hgc.cassetteTile(hid.iphi());
0180             edm::LogVerbatim(pid) << "Hit[" << all << ":" << allSc << "] " << hid << " Cassette " << cassette << " at ("
0181                                   << pos.x() << ", " << pos.y() << ", " << pos.z() << ") Valid " << valid1 << ":"
0182                                   << valid2;
0183           }
0184         } else {
0185           st1 << std::hex << id.rawId() << std::dec;
0186         }
0187         if (!valid) {
0188           edm::LogVerbatim("HGCalError") << "Invalid ID " << st1.str();
0189           ++bad;
0190         }
0191       }
0192     }
0193   }
0194   edm::LogVerbatim("HGCalSim") << "Total hits = " << all << ":" << nhits << " Good Silicon DetIds = " << allSi << ":"
0195                                << good << " Scintitllator = " << allSc << " Invalid = " << bad;
0196 }
0197 
0198 //define this as a plug-in
0199 DEFINE_FWK_MODULE(HGCalTestPartialWaferHits);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!