HGCalValidation/plugins/HGCGeometryValidation.cc

0001 #include <memory>
0002 #include <iostream>
0003
0004 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0005 #include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h"
0006
0007 #include "DataFormats/Common/interface/Handle.h"
0008 #include "DataFormats/ForwardDetId/interface/HGCSiliconDetId.h"
0009 #include "DataFormats/ForwardDetId/interface/HGCScintillatorDetId.h"
0010
0011 #include "FWCore/Framework/interface/Frameworkfwd.h"
0012 #include "FWCore/Framework/interface/Event.h"
0013 #include "FWCore/Framework/interface/MakerMacros.h"
0014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0015 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0016 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0017 #include "FWCore/ServiceRegistry/interface/Service.h"
0018 #include "FWCore/Utilities/interface/Exception.h"
0019 #include "FWCore/Utilities/interface/EDGetToken.h"
0020 #include "FWCore/Utilities/interface/transform.h"
0021 #include "DQMServices/Core/interface/DQMEDAnalyzer.h"
0022 #include "DQMServices/Core/interface/DQMStore.h"
0023
0024 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0025 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0026 #include "SimDataFormats/ValidationFormats/interface/PHGCalValidInfo.h"
0027
0028 #include "PhysicsTools/HepMCCandAlgos/interface/GenParticlesHelper.h"
0029
0030 const double mmtocm = 0.1;
0031
0032 class HGCGeometryValidation : public DQMEDAnalyzer {
0033 public:
0034   explicit HGCGeometryValidation(const edm::ParameterSet &);
0035   ~HGCGeometryValidation() override = default;
0036   static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0037
0038 protected:
0039   void dqmBeginRun(edm::Run const &, edm::EventSetup const &) override;
0040   void bookHistograms(DQMStore::IBooker &, edm::Run const &, edm::EventSetup const &) override;
0041   void analyze(const edm::Event &, const edm::EventSetup &) override;
0042
0043 private:
0044   const edm::EDGetTokenT<PHGCalValidInfo> g4Token_;
0045   const std::vector<std::string> geometrySource_;
0046   const int verbosity_;
0047   const std::vector<edm::ESGetToken<HGCalDDDConstants, IdealGeometryRecord>> geomToken_;
0048
0049   //HGCal geometry scheme
0050   std::vector<const HGCalDDDConstants *> hgcGeometry_;
0051
0052   //histogram related stuff
0053   MonitorElement *heedzVsZ, *heedyVsY, *heedxVsX;
0054   MonitorElement *hefdzVsZ, *hefdyVsY, *hefdxVsX;
0055   MonitorElement *hebdzVsZ, *hebdyVsY, *hebdxVsX;
0056   MonitorElement *heedzVsLayer, *hefdzVsLayer, *hebdzVsLayer;
0057   MonitorElement *heedyVsLayer, *hefdyVsLayer, *hebdyVsLayer;
0058   MonitorElement *heedxVsLayer, *hefdxVsLayer, *hebdxVsLayer;
0059   MonitorElement *heeXG4VsId, *hefXG4VsId, *hebXG4VsId;
0060   MonitorElement *heeYG4VsId, *hefYG4VsId, *hebYG4VsId;
0061   MonitorElement *heeZG4VsId, *hefZG4VsId, *hebZG4VsId;
0062   MonitorElement *hebLayerVsEnStep, *hefLayerVsEnStep, *heeLayerVsEnStep;
0063
0064   MonitorElement *heeTotEdepStep, *hefTotEdepStep, *hebTotEdepStep;
0065   MonitorElement *heedX, *heedY, *heedZ;
0066   MonitorElement *hefdX, *hefdY, *hefdZ;
0067   MonitorElement *hebdX, *hebdY, *hebdZ;
0068 };
0069
0070 HGCGeometryValidation::HGCGeometryValidation(const edm::ParameterSet &cfg)
0071     : g4Token_(consumes<PHGCalValidInfo>(cfg.getParameter<edm::InputTag>("g4Source"))),
0072       geometrySource_(cfg.getUntrackedParameter<std::vector<std::string>>("geometrySource")),
0073       verbosity_(cfg.getUntrackedParameter<int>("verbosity", 0)),
0074       geomToken_{edm::vector_transform(geometrySource_, [this](const std::string &name) {
0075         return esConsumes<HGCalDDDConstants, IdealGeometryRecord, edm::Transition::BeginRun>(edm::ESInputTag{"", name});
0076       })} {}
0077
0078 void HGCGeometryValidation::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0079   edm::ParameterSetDescription desc;
0080   desc.setUnknown();
0081   descriptions.addDefault(desc);
0082 }
0083
0084 void HGCGeometryValidation::dqmBeginRun(const edm::Run &iRun, const edm::EventSetup &iSetup) {
0085   for (size_t i = 0; i < geometrySource_.size(); i++) {
0086     hgcGeometry_.emplace_back(&iSetup.getData(geomToken_[i]));
0087   }
0088 }
0089
0090 void HGCGeometryValidation::bookHistograms(DQMStore::IBooker &iB, edm::Run const &, edm::EventSetup const &) {
0091   iB.setCurrentFolder("HGCAL/HGCalSimHitsV/Geometry");
0092
0093   //initiating histograms
0094   heeTotEdepStep = iB.book1D("heeTotEdepStep", "", 100, 0, 100);
0095   hefTotEdepStep = iB.book1D("hefTotEdepStep", "", 100, 0, 100);
0096   hebTotEdepStep = iB.book1D("hebTotEdepStep", "", 100, 0, 100);
0097
0098   hebLayerVsEnStep = iB.book2D("hebLayerVsEnStep", "", 25, 0, 25, 100, 0, 0.01);
0099   hefLayerVsEnStep = iB.book2D("hefLayerVsEnStep", "", 36, 0, 36, 100, 0, 0.01);
0100   heeLayerVsEnStep = iB.book2D("heeLayerVsEnStep", "", 84, 0, 84, 100, 0, 0.01);
0101
0102   heeXG4VsId = iB.book2D("heeXG4VsId", "", 600, -300, 300, 600, -300, 300);
0103   heeYG4VsId = iB.book2D("heeYG4VsId", "", 600, -300, 300, 600, -300, 300);
0104   heeZG4VsId = iB.book2D("heeZG4VsId", "", 3000, 320, 350, 3000, 320, 350);
0105
0106   hefXG4VsId = iB.book2D("hefXG4VsId", "", 600, -300, 300, 600, -300, 300);
0107   hefYG4VsId = iB.book2D("hefYG4VsId", "", 600, -300, 300, 600, -300, 300);
0108   hefZG4VsId = iB.book2D("hefZG4VsId", "", 6000, 350, 410, 6000, 350, 410);
0109
0110   hebXG4VsId = iB.book2D("hebXG4VsId", "", 600, -300, 300, 600, -300, 300);
0111   hebYG4VsId = iB.book2D("hebYG4VsId", "", 600, -300, 300, 600, -300, 300);
0112   hebZG4VsId = iB.book2D("hebZG4VsId", "", 220, 400, 620, 220, 400, 620);
0113
0114   heedzVsZ = iB.book2D("heedzVsZ", "", 600, 320, 350, 100, -1, 1);
0115   heedyVsY = iB.book2D("heedyVsY", "", 400, -200, 200, 100, -1, 1);
0116   heedxVsX = iB.book2D("heedxVsX", "", 400, -200, 200, 100, -1, 1);
0117
0118   hefdzVsZ = iB.book2D("hefdzVsZ", "", 1200, 350, 410, 100, -1, 1);
0119   hefdyVsY = iB.book2D("hefdyVsY", "", 400, -200, 200, 100, -1, 1);
0120   hefdxVsX = iB.book2D("hefdxVsX", "", 400, -200, 200, 100, -1, 1);
0121
0122   hebdzVsZ = iB.book2D("hebdzVsZ", "", 220, 400, 620, 100, -5, 5);
0123   hebdyVsY = iB.book2D("hebdyVsY", "", 400, -200, 200, 100, -5, 5);
0124   hebdxVsX = iB.book2D("hebdxVsX", "", 400, -200, 200, 100, -5, 5);
0125
0126   heedzVsLayer = iB.book2D("heedzVsLayer", "", 100, 0, 100, 100, -1, 1);
0127   hefdzVsLayer = iB.book2D("hefdzVsLayer", "", 40, 0, 40, 100, -1, 1);
0128   hebdzVsLayer = iB.book2D("hebdzVsLayer", "", 50, 0, 25, 100, -5, 5);
0129
0130   heedyVsLayer = iB.book2D("heedyVsLayer", "", 100, 0, 100, 100, -1, 1);
0131   hefdyVsLayer = iB.book2D("hefdyVsLayer", "", 40, 0, 40, 100, -1, 1);
0132   hebdyVsLayer = iB.book2D("hebdyVsLayer", "", 50, 0, 25, 100, -5, 5);
0133
0134   heedxVsLayer = iB.book2D("heedxVsLayer", "", 100, 0, 100, 100, -1, 1);
0135   hefdxVsLayer = iB.book2D("hefdxVsLayer", "", 40, 0, 40, 500, -1, 1);
0136   hebdxVsLayer = iB.book2D("hebdxVsLayer", "", 50, 0, 25, 500, -5, 5.0);
0137
0138   heedX = iB.book1D("heedX", "", 100, -1, 1);
0139   heedY = iB.book1D("heedY", "", 100, -1, 1);
0140   heedZ = iB.book1D("heedZ", "", 100, -1, 1);
0141
0142   hefdX = iB.book1D("hefdX", "", 100, -1, 1);
0143   hefdY = iB.book1D("hefdY", "", 100, -1, 1);
0144   hefdZ = iB.book1D("hefdZ", "", 100, -1, 1);
0145
0146   hebdX = iB.book1D("hebdX", "", 100, -1, 1);
0147   hebdY = iB.book1D("hebdY", "", 100, -1, 1);
0148   hebdZ = iB.book1D("hebdZ", "", 100, -1, 1);
0149 }
0150
0151 void HGCGeometryValidation::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0152   //Accessing G4 information
0153   const edm::Handle<PHGCalValidInfo> &infoLayer = iEvent.getHandle(g4Token_);
0154
0155   if (infoLayer.isValid()) {
0156     //step vertex information
0157     std::vector<float> hitVtxX = infoLayer->hitvtxX();
0158     std::vector<float> hitVtxY = infoLayer->hitvtxY();
0159     std::vector<float> hitVtxZ = infoLayer->hitvtxZ();
0160     const std::vector<unsigned int> &hitDet = infoLayer->hitDets();
0161     const std::vector<unsigned int> &hitIdx = infoLayer->hitIndex();
0162
0163     //energy information
0164     const std::vector<float> &edepLayerEE = infoLayer->eehgcEdep();
0165     const std::vector<float> &edepLayerHE = infoLayer->hefhgcEdep();
0166     const std::vector<float> &edepLayerHB = infoLayer->hebhgcEdep();
0167
0168     unsigned int i;
0169     for (i = 0; i < edepLayerEE.size(); i++) {
0170       heeLayerVsEnStep->Fill(i, edepLayerEE[i]);
0171     }
0172
0173     for (i = 0; i < edepLayerHE.size(); i++) {
0174       hefLayerVsEnStep->Fill(i, edepLayerHE[i]);
0175     }
0176
0177     for (i = 0; i < edepLayerHB.size(); i++) {
0178       hebLayerVsEnStep->Fill(i, edepLayerHB[i]);
0179     }
0180
0181     //fill total energy deposited
0182     heeTotEdepStep->Fill((double)infoLayer->eeTotEdep());
0183     hefTotEdepStep->Fill((double)infoLayer->hefTotEdep());
0184     hebTotEdepStep->Fill((double)infoLayer->hebTotEdep());
0185
0186     //loop over all hits
0187     for (unsigned int i = 0; i < hitVtxX.size(); i++) {
0188       hitVtxX[i] *= mmtocm;
0189       hitVtxY[i] *= mmtocm;
0190       hitVtxZ[i] *= mmtocm;
0191
0192       double xx, yy;
0193       int dtype(0), layer(0), zside(1);
0194       std::pair<float, float> xy;
0195       if ((hitDet[i] == static_cast<unsigned int>(DetId::HGCalEE)) ||
0196           (hitDet[i] == static_cast<unsigned int>(DetId::HGCalHSi))) {
0197         HGCSiliconDetId id(hitIdx[i]);
0198         dtype = (id.det() == DetId::HGCalEE) ? 0 : 1;
0199         layer = id.layer();
0200         zside = id.zside();
0201         xy = hgcGeometry_[dtype]->locateCell(
0202             zside, layer, id.waferU(), id.waferV(), id.cellU(), id.cellV(), true, true, false, false, false);
0203       } else {
0204         HGCScintillatorDetId id(hitIdx[i]);
0205         dtype = 2;
0206         layer = id.layer();
0207         zside = id.zside();
0208         xy = hgcGeometry_[dtype]->locateCellTrap(zside, layer, id.ietaAbs(), id.iphi(), true, false);
0209       }
0210       double zz = hgcGeometry_[dtype]->waferZ(layer, true);  //cm
0211       if (zside < 0)
0212         zz = -zz;
0213       xx = (zside < 0) ? -xy.first : xy.first;
0214       yy = xy.second;
0215
0216       if (dtype == 0) {
0217         heedzVsZ->Fill(zz, (hitVtxZ[i] - zz));
0218         heedyVsY->Fill(yy, (hitVtxY[i] - yy));
0219         heedxVsX->Fill(xx, (hitVtxX[i] - xx));
0220
0221         heeXG4VsId->Fill(hitVtxX[i], xx);
0222         heeYG4VsId->Fill(hitVtxY[i], yy);
0223         heeZG4VsId->Fill(hitVtxZ[i], zz);
0224
0225         heedzVsLayer->Fill(layer, (hitVtxZ[i] - zz));
0226         heedyVsLayer->Fill(layer, (hitVtxY[i] - yy));
0227         heedxVsLayer->Fill(layer, (hitVtxX[i] - xx));
0228
0229         heedX->Fill((hitVtxX[i] - xx));
0230         heedZ->Fill((hitVtxZ[i] - zz));
0231         heedY->Fill((hitVtxY[i] - yy));
0232
0233       } else if (dtype == 1) {
0234         hefdzVsZ->Fill(zz, (hitVtxZ[i] - zz));
0235         hefdyVsY->Fill(yy, (hitVtxY[i] - yy));
0236         hefdxVsX->Fill(xx, (hitVtxX[i] - xx));
0237
0238         hefXG4VsId->Fill(hitVtxX[i], xx);
0239         hefYG4VsId->Fill(hitVtxY[i], yy);
0240         hefZG4VsId->Fill(hitVtxZ[i], zz);
0241
0242         hefdzVsLayer->Fill(layer, (hitVtxZ[i] - zz));
0243         hefdyVsLayer->Fill(layer, (hitVtxY[i] - yy));
0244         hefdxVsLayer->Fill(layer, (hitVtxX[i] - xx));
0245
0246         hefdX->Fill((hitVtxX[i] - xx));
0247         hefdZ->Fill((hitVtxZ[i] - zz));
0248         hefdY->Fill((hitVtxY[i] - yy));
0249
0250       } else {
0251         hebdzVsZ->Fill(zz, (hitVtxZ[i] - zz));
0252         hebdyVsY->Fill(yy, (hitVtxY[i] - yy));
0253         hebdxVsX->Fill(xx, (hitVtxX[i] - xx));
0254
0255         hebXG4VsId->Fill(hitVtxX[i], xx);
0256         hebYG4VsId->Fill(hitVtxY[i], yy);
0257         hebZG4VsId->Fill(hitVtxZ[i], zz);
0258
0259         hebdzVsLayer->Fill(layer, (hitVtxZ[i] - zz));
0260         hebdyVsLayer->Fill(layer, (hitVtxY[i] - yy));
0261         hebdxVsLayer->Fill(layer, (hitVtxX[i] - xx));
0262
0263         hebdX->Fill((hitVtxX[i] - xx));
0264         hebdZ->Fill((hitVtxZ[i] - zz));
0265         hebdY->Fill((hitVtxY[i] - yy));
0266       }
0267     }  //end G4 hits
0268
0269   } else {
0270     if (verbosity_ > 0)
0271       edm::LogVerbatim("HGCalValid") << "HGCGeometryValidation::No PHGCalInfo";
0272   }
0273 }
0274
0275 //define this as a plug-in
0276 DEFINE_FWK_MODULE(HGCGeometryValidation);