HGCalValidation/test/HGCGeometryCheck.cc

0001 #include <memory>
0002 #include <iostream>
0003
0004 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0005
0006 #include "DataFormats/ForwardDetId/interface/HGCSiliconDetId.h"
0007 #include "DataFormats/ForwardDetId/interface/HGCScintillatorDetId.h"
0008
0009 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0010 #include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h"
0011
0012 #include "FWCore/Framework/interface/Frameworkfwd.h"
0013 #include "FWCore/Framework/interface/Event.h"
0014 #include "FWCore/Framework/interface/MakerMacros.h"
0015 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0016 #include "FWCore/Framework/interface/ESHandle.h"
0017 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0018 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0019 #include "FWCore/ServiceRegistry/interface/Service.h"
0020 #include "FWCore/Utilities/interface/EDGetToken.h"
0021 #include "FWCore/Utilities/interface/transform.h"
0022
0023 #include "SimDataFormats/ValidationFormats/interface/PHGCalValidInfo.h"
0024 #include "SimDataFormats/CaloTest/interface/HGCalTestNumbering.h"
0025
0026 #include <TH2.h>
0027
0028 //#define EDM_ML_DEBUG
0029
0030 class HGCGeometryCheck : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0031 public:
0032   explicit HGCGeometryCheck(const edm::ParameterSet &);
0033   ~HGCGeometryCheck() override = default;
0034   static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0035
0036 protected:
0037   void beginJob() override;
0038   void beginRun(edm::Run const &, edm::EventSetup const &) override;
0039   void analyze(edm::Event const &, edm::EventSetup const &) override;
0040   void endRun(edm::Run const &, edm::EventSetup const &) override {}
0041   virtual void beginLuminosityBlock(edm::LuminosityBlock const &, edm::EventSetup const &) {}
0042   virtual void endLuminosityBlock(edm::LuminosityBlock const &, edm::EventSetup const &) {}
0043
0044 private:
0045   const edm::EDGetTokenT<PHGCalValidInfo> g4Token_;
0046   const std::vector<std::string> geometrySource_;
0047   const std::vector<edm::ESGetToken<HGCalDDDConstants, IdealGeometryRecord> > tok_hgcGeom_;
0048
0049   //HGCal geometry scheme
0050   std::vector<const HGCalDDDConstants *> hgcGeometry_;
0051
0052   //histogram related stuff
0053   TH2F *heedzVsZ, *hefdzVsZ, *hebdzVsZ;
0054   TH2F *heezVsLayer, *hefzVsLayer, *hebzVsLayer;
0055   TH2F *heerVsLayer, *hefrVsLayer, *hebrVsLayer;
0056
0057   static constexpr double mmTocm_ = 0.1;
0058 };
0059
0060 HGCGeometryCheck::HGCGeometryCheck(const edm::ParameterSet &cfg)
0061     : g4Token_(consumes<PHGCalValidInfo>(cfg.getParameter<edm::InputTag>("g4Source"))),
0062       geometrySource_(cfg.getUntrackedParameter<std::vector<std::string> >("geometrySource")),
0063       tok_hgcGeom_{edm::vector_transform(geometrySource_, [this](const std::string &name) {
0064         return esConsumes<HGCalDDDConstants, IdealGeometryRecord, edm::Transition::BeginRun>(edm::ESInputTag{"", name});
0065       })} {
0066   usesResource(TFileService::kSharedResource);
0067
0068   edm::LogVerbatim("HGCalValid") << "HGCGeometryCheck:: use information from "
0069                                  << cfg.getParameter<edm::InputTag>("g4Source") << " and " << geometrySource_.size()
0070                                  << " geometry records:";
0071   for (unsigned int k = 0; k < geometrySource_.size(); ++k)
0072     edm::LogVerbatim("HGCalValid") << "[ " << k << "] " << geometrySource_[k];
0073 }
0074
0075 void HGCGeometryCheck::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0076   std::vector<std::string> names = {"HGCalEESensitive", "HGCalHESiliconSensitive", "HGCalHEScintillatorSensitive"};
0077   edm::ParameterSetDescription desc;
0078   desc.addUntracked<std::vector<std::string> >("geometrySource", names);
0079   desc.add<edm::InputTag>("g4Source", edm::InputTag("g4SimHits", "HGCalInfoLayer"));
0080   descriptions.add("hgcGeomCheck", desc);
0081 }
0082
0083 void HGCGeometryCheck::beginJob() {
0084   //initiating fileservice
0085   edm::Service<TFileService> fs;
0086
0087   //initiating histograms
0088   heedzVsZ = fs->make<TH2F>("heedzVsZ", "", 1400, 315, 385, 100, -1, 1);
0089   hefdzVsZ = fs->make<TH2F>("hefdzVsZ", "", 2000, 350, 550, 100, -1, 1);
0090   hebdzVsZ = fs->make<TH2F>("hebdzVsZ", "", 360, 380, 560, 100, -5, 5);
0091
0092   heezVsLayer = fs->make<TH2F>("heezVsLayer", "", 100, 0, 100, 1400, 315, 385);
0093   hefzVsLayer = fs->make<TH2F>("hefzVsLayer", "", 40, 0, 40, 2000, 350, 550);
0094   hebzVsLayer = fs->make<TH2F>("hebzVsLayer", "", 50, 0, 25, 360, 380, 560);
0095
0096   heerVsLayer = fs->make<TH2F>("heerVsLayer", "", 100, 0, 100, 600, 0, 300);
0097   hefrVsLayer = fs->make<TH2F>("hefrVsLayer", "", 40, 0, 40, 600, 0, 300);
0098   hebrVsLayer = fs->make<TH2F>("hebrVsLayer", "", 50, 0, 25, 600, 0, 300);
0099 }
0100
0101 void HGCGeometryCheck::beginRun(const edm::Run &, const edm::EventSetup &iSetup) {
0102   //initiating hgc geometry
0103   for (size_t i = 0; i < geometrySource_.size(); i++) {
0104     const edm::ESHandle<HGCalDDDConstants> &hgcGeom = iSetup.getHandle(tok_hgcGeom_[i]);
0105     if (hgcGeom.isValid()) {
0106       hgcGeometry_.push_back(hgcGeom.product());
0107       edm::LogVerbatim("HGCalValid") << "Initialize geometry for " << geometrySource_[i];
0108     } else {
0109       edm::LogWarning("HGCalValid") << "Cannot initiate HGCalDDDConstants for " << geometrySource_[i];
0110     }
0111   }
0112 }
0113
0114 void HGCGeometryCheck::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0115 #ifdef EDM_ML_DEBUG
0116   edm::LogVerbatim("HGCalValid") << "HGCGeometryCheck::Run " << iEvent.id().run() << " Event " << iEvent.id().event()
0117                                  << " Luminosity " << iEvent.luminosityBlock() << " Bunch " << iEvent.bunchCrossing();
0118 #endif
0119   //Accessing G4 information
0120   const edm::Handle<PHGCalValidInfo> &infoLayer = iEvent.getHandle(g4Token_);
0121
0122   if (infoLayer.isValid()) {
0123     //step vertex information
0124     const std::vector<float> &hitVtxX = infoLayer->hitvtxX();
0125     const std::vector<float> &hitVtxY = infoLayer->hitvtxY();
0126     const std::vector<float> &hitVtxZ = infoLayer->hitvtxZ();
0127     const std::vector<unsigned int> &hitDet = infoLayer->hitDets();
0128     const std::vector<unsigned int> &hitIdx = infoLayer->hitIndex();
0129
0130     //loop over all hits
0131     for (unsigned int i = 0; i < hitVtxZ.size(); i++) {
0132       double xx = mmTocm_ * hitVtxX[i];
0133       double yy = mmTocm_ * hitVtxY[i];
0134       double zz = mmTocm_ * hitVtxZ[i];
0135       double rr = sqrt(xx * xx + yy * yy);
0136       if ((hitDet[i] == static_cast<unsigned int>(DetId::Forward)) ||
0137           (hitDet[i] == static_cast<unsigned int>(DetId::HGCalEE)) ||
0138           (hitDet[i] == static_cast<unsigned int>(DetId::HGCalHSi)) ||
0139           (hitDet[i] == static_cast<unsigned int>(DetId::HGCalHSc))) {
0140         int dtype(0), layer(0), zside(1);
0141         if ((hitDet[i] == static_cast<unsigned int>(DetId::HGCalEE)) ||
0142             (hitDet[i] == static_cast<unsigned int>(DetId::HGCalHSi))) {
0143           HGCSiliconDetId id(hitIdx[i]);
0144           dtype = (id.det() == DetId::HGCalEE) ? 0 : 1;
0145           layer = id.layer();
0146           zside = id.zside();
0147         } else {
0148           HGCScintillatorDetId id(hitIdx[i]);
0149           dtype = 2;
0150           layer = id.layer();
0151           zside = id.zside();
0152         }
0153         double zp = hgcGeometry_[dtype]->waferZ(layer, true);  //cm
0154         if (zside < 0)
0155           zp = -zp;
0156 #ifdef EDM_ML_DEBUG
0157         edm::LogVerbatim("HGCalValid") << "Info[" << i << "] Detector Information " << hitDet[i] << ":" << zside << ":"
0158                                        << layer << " Z " << zp << ":" << zz << " R " << rr;
0159 #endif
0160         if (dtype == 0) {
0161           heedzVsZ->Fill(zp, (zz - zp));
0162           heezVsLayer->Fill(layer, zz);
0163           heerVsLayer->Fill(layer, rr);
0164         } else if (dtype == 1) {
0165           hefdzVsZ->Fill(zp, (zz - zp));
0166           hefzVsLayer->Fill(layer, zz);
0167           hefrVsLayer->Fill(layer, rr);
0168         } else {
0169           hebdzVsZ->Fill(zp, (zz - zp));
0170           hebzVsLayer->Fill(layer, zz);
0171           hebrVsLayer->Fill(layer, rr);
0172         }
0173       }
0174     }  //end G4 hits
0175   } else {
0176     edm::LogWarning("HGCalValid") << "No PHGCalInfo " << std::endl;
0177   }
0178 }
0179
0180 //define this as a plug-in
0181 DEFINE_FWK_MODULE(HGCGeometryCheck);