Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-09-22 23:03:58

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