Calo/plugins/HGCalTestScintHits.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/HGCScintillatorDetId.h"
0016
0017 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0018 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0019 #include "SimG4CMS/Calo/interface/CaloSimUtils.h"
0020
0021 #include "Geometry/CaloTopology/interface/HGCalTopology.h"
0022 #include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
0023 #include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h"
0024 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0025
0026 #include <fstream>
0027 #include <string>
0028 #include <vector>
0029
0030 class HGCalTestScintHits : public edm::one::EDAnalyzer<> {
0031 public:
0032   HGCalTestScintHits(const edm::ParameterSet& ps);
0033   ~HGCalTestScintHits() override = default;
0034   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0035
0036 protected:
0037   void analyze(edm::Event const&, edm::EventSetup const&) override;
0038   void beginJob() override {}
0039   void endJob() override {}
0040
0041 private:
0042   const std::string g4Label_, caloHitSource_, nameSense_, tileFileName_;
0043   const edm::EDGetTokenT<edm::PCaloHitContainer> tok_calo_;
0044   const edm::ESGetToken<HGCalGeometry, IdealGeometryRecord> geomToken_;
0045   std::vector<int> tiles_;
0046 };
0047
0048 HGCalTestScintHits::HGCalTestScintHits(const edm::ParameterSet& ps)
0049     : g4Label_(ps.getParameter<std::string>("moduleLabel")),
0050       caloHitSource_(ps.getParameter<std::string>("caloHitSource")),
0051       nameSense_(ps.getParameter<std::string>("nameSense")),
0052       tileFileName_(ps.getParameter<std::string>("tileFileName")),
0053       tok_calo_(consumes<edm::PCaloHitContainer>(edm::InputTag(g4Label_, caloHitSource_))),
0054       geomToken_(esConsumes<HGCalGeometry, IdealGeometryRecord>(edm::ESInputTag{"", nameSense_})) {
0055   edm::LogVerbatim("HGCalSim") << "Test Hit ID using SimHits for " << nameSense_ << " with module Label: " << g4Label_
0056                                << "   Hits: " << caloHitSource_ << " Tile file " << tileFileName_;
0057   if (!tileFileName_.empty()) {
0058     edm::FileInPath filetmp("SimG4CMS/Calo/data/" + tileFileName_);
0059     std::string fileName = filetmp.fullPath();
0060     std::ifstream fInput(fileName.c_str());
0061     if (!fInput.good()) {
0062       edm::LogVerbatim("HGCalSim") << "Cannot open file " << fileName;
0063     } else {
0064       char buffer[80];
0065       while (fInput.getline(buffer, 80)) {
0066         std::vector<std::string> items = CaloSimUtils::splitString(std::string(buffer));
0067         if (items.size() > 2) {
0068           int layer = std::atoi(items[0].c_str());
0069           int ring = std::atoi(items[1].c_str());
0070           int phi = std::atoi(items[2].c_str());
0071           tiles_.emplace_back(HGCalTileIndex::tileIndex(layer, ring, phi));
0072         }
0073       }
0074       edm::LogVerbatim("HGCalSim") << "Reads in " << tiles_.size() << " tile information from " << fileName;
0075       fInput.close();
0076     }
0077   }
0078 }
0079
0080 void HGCalTestScintHits::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0081   edm::ParameterSetDescription desc;
0082   desc.add<std::string>("moduleLabel", "g4SimHits");
0083   desc.add<std::string>("caloHitSource", "HGCHitsHEback");
0084   desc.add<std::string>("nameSense", "HGCalHEScintillatorSensitive");
0085   desc.add<std::string>("tileFileName", "");
0086   descriptions.add("hgcalHitScintillator", desc);
0087 }
0088
0089 void HGCalTestScintHits::analyze(const edm::Event& e, const edm::EventSetup& iS) {
0090   // get hcalGeometry
0091   const HGCalGeometry* geom = &iS.getData(geomToken_);
0092   const HGCalDDDConstants& hgc = geom->topology().dddConstants();
0093   int firstLayer = hgc.firstLayer() - 1;
0094   // get the hit collection
0095   const edm::Handle<edm::PCaloHitContainer>& hitsCalo = e.getHandle(tok_calo_);
0096   bool getHits = (hitsCalo.isValid());
0097   uint32_t nhits = (getHits) ? hitsCalo->size() : 0;
0098   uint32_t good(0), all(0);
0099   edm::LogVerbatim("HGCalSim") << "HGCalTestScintHits: Input flags Hits " << getHits << " with " << nhits
0100                                << " hits first Layer " << firstLayer;
0101
0102   if (getHits) {
0103     std::vector<PCaloHit> hits;
0104     hits.insert(hits.end(), hitsCalo->begin(), hitsCalo->end());
0105     if (!hits.empty()) {
0106       // Loop over all hits
0107       for (auto hit : hits) {
0108         ++all;
0109         DetId id(hit.id());
0110         HGCScintillatorDetId hid(id);
0111         std::pair<int, int> typm = hgc.tileType(hid.layer(), hid.ring(), 0);
0112         if (typm.first >= 0) {
0113           hid.setType(typm.first);
0114           hid.setSiPM(typm.second);
0115           id = static_cast<DetId>(hid);
0116         }
0117         bool toCheck(true);
0118         if (!tiles_.empty()) {
0119           int indx = HGCalTileIndex::tileIndex(firstLayer + hid.layer(), hid.ring(), hid.iphi());
0120           if (std::find(tiles_.begin(), tiles_.end(), indx) != tiles_.end())
0121             toCheck = true;
0122         }
0123         if (toCheck) {
0124           ++good;
0125           GlobalPoint pos = geom->getPosition(id);
0126           bool valid1 = geom->topology().valid(id);
0127           bool valid2 = hgc.isValidTrap(hid.zside(), hid.layer(), hid.ring(), hid.iphi());
0128           edm::LogVerbatim("HGCalSim") << "Hit[" << all << ":" << good << "]" << hid << " at (" << pos.x() << ", "
0129                                        << pos.y() << ", " << pos.z() << ") Validity " << valid1 << ":" << valid2
0130                                        << " Energy " << hit.energy() << " Time " << hit.time();
0131         }
0132       }
0133     }
0134   }
0135   edm::LogVerbatim("HGCalSim") << "Total hits = " << all << ":" << nhits << " Good DetIds = " << good;
0136 }
0137
0138 //define this as a plug-in
0139 DEFINE_FWK_MODULE(HGCalTestScintHits);