Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-10-10 04:02:15

0001 //// -*- C++ -*-
0002 //
0003 // Package:    HGCalHitValidation
0004 // Class:      HGCalHitValidation
0005 //
0006 /**\class HGCalHitValidation HGCalHitValidation.cc Validation/HGCalValidation/plugins/HGCalHitValidation.cc
0007 
0008  Description: [one line class summary]
0009 
0010  Implementation:
0011     [Notes on implementation]
0012 */
0013 
0014 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0015 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0016 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0017 #include "Geometry/HGCalCommonData/interface/HGCalGeometryMode.h"
0018 #include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
0019 #include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h"
0020 
0021 #include "DataFormats/Common/interface/Handle.h"
0022 #include "DataFormats/HGCRecHit/interface/HGCRecHit.h"
0023 #include "DataFormats/HGCRecHit/interface/HGCRecHitCollections.h"
0024 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0025 
0026 #include "FWCore/Framework/interface/Frameworkfwd.h"
0027 #include "FWCore/Framework/interface/Event.h"
0028 #include "FWCore/Framework/interface/ESHandle.h"
0029 #include "FWCore/ServiceRegistry/interface/Service.h"
0030 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0031 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0032 #include "FWCore/Utilities/interface/Exception.h"
0033 #include "FWCore/Utilities/interface/EDGetToken.h"
0034 #include "FWCore/Utilities/interface/transform.h"
0035 #include "DQMServices/Core/interface/DQMEDAnalyzer.h"
0036 #include "DQMServices/Core/interface/DQMStore.h"
0037 
0038 #include "SimG4CMS/Calo/interface/HGCNumberingScheme.h"
0039 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0040 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0041 #include "SimDataFormats/CaloTest/interface/HGCalTestNumbering.h"
0042 
0043 #include <cmath>
0044 #include <memory>
0045 #include <iostream>
0046 #include <string>
0047 #include <vector>
0048 
0049 //#define EDM_ML_DEBUG
0050 
0051 class HGCalHitValidation : public DQMEDAnalyzer {
0052 public:
0053   explicit HGCalHitValidation(const edm::ParameterSet&);
0054   ~HGCalHitValidation() override = default;
0055   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0056 
0057 protected:
0058   typedef std::tuple<float, GlobalPoint> HGCHitTuple;
0059 
0060   void dqmBeginRun(edm::Run const&, edm::EventSetup const&) override;
0061   void bookHistograms(DQMStore::IBooker&, edm::Run const&, edm::EventSetup const&) override;
0062   void analyze(const edm::Event&, const edm::EventSetup&) override;
0063   void analyzeHGCalSimHit(edm::Handle<std::vector<PCaloHit>> const& simHits,
0064                           int idet,
0065                           MonitorElement* hist,
0066                           std::map<unsigned int, HGCHitTuple>&);
0067 
0068 private:
0069   //HGC Geometry
0070   std::vector<const HGCalDDDConstants*> hgcCons_;
0071   std::vector<const HGCalGeometry*> hgcGeometry_;
0072   const std::vector<std::string> geometrySource_;
0073   const std::vector<int> ietaExcludeBH_;
0074   const edm::EDGetTokenT<std::vector<PCaloHit>> eeSimHitToken_;
0075   const edm::EDGetTokenT<std::vector<PCaloHit>> fhSimHitToken_;
0076   const edm::EDGetTokenT<std::vector<PCaloHit>> bhSimHitToken_;
0077   const edm::EDGetTokenT<HGCeeRecHitCollection> eeRecHitToken_;
0078   const edm::EDGetTokenT<HGChefRecHitCollection> fhRecHitToken_;
0079   const edm::EDGetTokenT<HGChebRecHitCollection> bhRecHitToken_;
0080   const std::vector<edm::ESGetToken<HGCalDDDConstants, IdealGeometryRecord>> tok_ddd_;
0081   const std::vector<edm::ESGetToken<HGCalGeometry, IdealGeometryRecord>> tok_geom_;
0082 
0083   //histogram related stuff
0084   MonitorElement *heedzVsZ, *heedyVsY, *heedxVsX;
0085   MonitorElement *hefdzVsZ, *hefdyVsY, *hefdxVsX;
0086   MonitorElement *hebdzVsZ, *hebdPhiVsPhi, *hebdEtaVsEta;
0087 
0088   MonitorElement *heeRecVsSimZ, *heeRecVsSimY, *heeRecVsSimX;
0089   MonitorElement *hefRecVsSimZ, *hefRecVsSimY, *hefRecVsSimX;
0090   MonitorElement *hebRecVsSimZ, *hebRecVsSimY, *hebRecVsSimX;
0091   MonitorElement *heeEnSimRec, *hefEnSimRec, *hebEnSimRec;
0092 
0093   MonitorElement *hebEnRec, *hebEnSim;
0094   MonitorElement *hefEnRec, *hefEnSim;
0095   MonitorElement *heeEnRec, *heeEnSim;
0096 };
0097 
0098 HGCalHitValidation::HGCalHitValidation(const edm::ParameterSet& cfg)
0099     : geometrySource_(cfg.getParameter<std::vector<std::string>>("geometrySource")),
0100       ietaExcludeBH_(cfg.getParameter<std::vector<int>>("ietaExcludeBH")),
0101       eeSimHitToken_(consumes<std::vector<PCaloHit>>(cfg.getParameter<edm::InputTag>("eeSimHitSource"))),
0102       fhSimHitToken_(consumes<std::vector<PCaloHit>>(cfg.getParameter<edm::InputTag>("fhSimHitSource"))),
0103       bhSimHitToken_(consumes<std::vector<PCaloHit>>(cfg.getParameter<edm::InputTag>("bhSimHitSource"))),
0104       eeRecHitToken_(consumes<HGCeeRecHitCollection>(cfg.getParameter<edm::InputTag>("eeRecHitSource"))),
0105       fhRecHitToken_(consumes<HGChefRecHitCollection>(cfg.getParameter<edm::InputTag>("fhRecHitSource"))),
0106       bhRecHitToken_(consumes<HGChebRecHitCollection>(cfg.getParameter<edm::InputTag>("bhRecHitSource"))),
0107       tok_ddd_{
0108           edm::vector_transform(geometrySource_,
0109                                 [this](const std::string& name) {
0110                                   return esConsumes<HGCalDDDConstants, IdealGeometryRecord, edm::Transition::BeginRun>(
0111                                       edm::ESInputTag{"", name});
0112                                 })},
0113       tok_geom_{edm::vector_transform(geometrySource_, [this](const std::string& name) {
0114         return esConsumes<HGCalGeometry, IdealGeometryRecord, edm::Transition::BeginRun>(edm::ESInputTag{"", name});
0115       })} {
0116 #ifdef EDM_ML_DEBUG
0117   edm::LogInfo("HGCalValid") << "Exclude the following " << ietaExcludeBH_.size() << " ieta values from BH plots";
0118   for (unsigned int k = 0; k < ietaExcludeBH_.size(); ++k)
0119     edm::LogInfo("HGCalValid") << " [" << k << "] " << ietaExcludeBH_[k];
0120 #endif
0121 }
0122 
0123 void HGCalHitValidation::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0124   edm::ParameterSetDescription desc;
0125   std::vector<std::string> source = {"HGCalEESensitive", "HGCalHESiliconSensitive", "HGCalHEScintillatorSensitive"};
0126   desc.add<std::vector<std::string>>("geometrySource", source);
0127   desc.add<edm::InputTag>("eeSimHitSource", edm::InputTag("g4SimHits", "HGCHitsEE"));
0128   desc.add<edm::InputTag>("fhSimHitSource", edm::InputTag("g4SimHits", "HGCHitsHEfront"));
0129   desc.add<edm::InputTag>("bhSimHitSource", edm::InputTag("g4SimHits", "HGCHitsHEback"));
0130   desc.add<edm::InputTag>("eeRecHitSource", edm::InputTag("HGCalRecHit", "HGCEERecHits"));
0131   desc.add<edm::InputTag>("fhRecHitSource", edm::InputTag("HGCalRecHit", "HGCHEFRecHits"));
0132   desc.add<edm::InputTag>("bhRecHitSource", edm::InputTag("HGCalRecHit", "HGCHEBRecHits"));
0133   std::vector<int> dummy;
0134   desc.add<std::vector<int>>("ietaExcludeBH", dummy);
0135   descriptions.add("hgcalHitValidation", desc);
0136 }
0137 
0138 void HGCalHitValidation::bookHistograms(DQMStore::IBooker& iB, edm::Run const&, edm::EventSetup const&) {
0139   iB.setCurrentFolder("HGCAL/HGCalSimHitsV/HitValidation");
0140 
0141   //initiating histograms
0142   heedzVsZ = iB.book2D("heedzVsZ", "", 720, -360, 360, 100, -0.1, 0.1);
0143   heedyVsY = iB.book2D("heedyVsY", "", 400, -200, 200, 100, -0.02, 0.02);
0144   heedxVsX = iB.book2D("heedxVsX", "", 400, -200, 200, 100, -0.02, 0.02);
0145   heeRecVsSimZ = iB.book2D("heeRecVsSimZ", "", 720, -360, 360, 720, -360, 360);
0146   heeRecVsSimY = iB.book2D("heeRecVsSimY", "", 400, -200, 200, 400, -200, 200);
0147   heeRecVsSimX = iB.book2D("heeRecVsSimX", "", 400, -200, 200, 400, -200, 200);
0148 
0149   hefdzVsZ = iB.book2D("hefdzVsZ", "", 820, -410, 410, 100, -0.1, 0.1);
0150   hefdyVsY = iB.book2D("hefdyVsY", "", 400, -200, 200, 100, -0.02, 0.02);
0151   hefdxVsX = iB.book2D("hefdxVsX", "", 400, -200, 200, 100, -0.02, 0.02);
0152   hefRecVsSimZ = iB.book2D("hefRecVsSimZ", "", 820, -410, 410, 820, -410, 410);
0153   hefRecVsSimY = iB.book2D("hefRecVsSimY", "", 400, -200, 200, 400, -200, 200);
0154   hefRecVsSimX = iB.book2D("hefRecVsSimX", "", 400, -200, 200, 400, -200, 200);
0155 
0156   hebdzVsZ = iB.book2D("hebdzVsZ", "", 1080, -540, 540, 100, -1.0, 1.0);
0157   hebdPhiVsPhi = iB.book2D("hebdPhiVsPhi", "", M_PI * 100, -0.5, M_PI + 0.5, 200, -0.2, 0.2);
0158   hebdEtaVsEta = iB.book2D("hebdEtaVsEta", "", 1000, -5, 5, 200, -0.1, 0.1);
0159   hebRecVsSimZ = iB.book2D("hebRecVsSimZ", "", 1080, -540, 540, 1080, -540, 540);
0160   hebRecVsSimY = iB.book2D("hebRecVsSimY", "", 400, -200, 200, 400, -200, 200);
0161   hebRecVsSimX = iB.book2D("hebRecVsSimX", "", 400, -200, 200, 400, -200, 200);
0162 
0163   heeEnRec = iB.book1D("heeEnRec", "", 1000, 0, 10);
0164   heeEnSim = iB.book1D("heeEnSim", "", 1000, 0, 0.01);
0165   heeEnSimRec = iB.book2D("heeEnSimRec", "", 200, 0, 0.002, 200, 0, 0.2);
0166 
0167   hefEnRec = iB.book1D("hefEnRec", "", 1000, 0, 10);
0168   hefEnSim = iB.book1D("hefEnSim", "", 1000, 0, 0.01);
0169   hefEnSimRec = iB.book2D("hefEnSimRec", "", 200, 0, 0.001, 200, 0, 0.5);
0170 
0171   hebEnRec = iB.book1D("hebEnRec", "", 1000, 0, 15);
0172   hebEnSim = iB.book1D("hebEnSim", "", 1000, 0, 0.01);
0173   hebEnSimRec = iB.book2D("hebEnSimRec", "", 200, 0, 0.02, 200, 0, 4);
0174 }
0175 
0176 void HGCalHitValidation::dqmBeginRun(edm::Run const& iRun, edm::EventSetup const& iSetup) {
0177   //initiating hgc Geometry
0178   for (size_t i = 0; i < geometrySource_.size(); i++) {
0179     const edm::ESHandle<HGCalDDDConstants>& hgcCons = iSetup.getHandle(tok_ddd_[i]);
0180     if (hgcCons.isValid()) {
0181       hgcCons_.push_back(hgcCons.product());
0182     } else {
0183       edm::LogWarning("HGCalValid") << "Cannot initiate HGCalDDDConstants for " << geometrySource_[i] << std::endl;
0184     }
0185     const edm::ESHandle<HGCalGeometry>& hgcGeom = iSetup.getHandle(tok_geom_[i]);
0186     if (hgcGeom.isValid()) {
0187       hgcGeometry_.push_back(hgcGeom.product());
0188     } else {
0189       edm::LogWarning("HGCalValid") << "Cannot initiate HGCalGeometry for " << geometrySource_[i] << std::endl;
0190     }
0191   }
0192 }
0193 
0194 void HGCalHitValidation::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0195   std::map<unsigned int, HGCHitTuple> eeHitRefs, fhHitRefs, bhHitRefs;
0196 
0197   //Accesing ee simhits
0198   const edm::Handle<std::vector<PCaloHit>>& eeSimHits = iEvent.getHandle(eeSimHitToken_);
0199 
0200   if (eeSimHits.isValid()) {
0201     analyzeHGCalSimHit(eeSimHits, 0, heeEnSim, eeHitRefs);
0202 #ifdef EDM_ML_DEBUG
0203     for (std::map<unsigned int, HGCHitTuple>::iterator itr = eeHitRefs.begin(); itr != eeHitRefs.end(); ++itr) {
0204       int idx = std::distance(eeHitRefs.begin(), itr);
0205       edm::LogInfo("HGCalValid") << "EEHit[" << idx << "] " << std::hex << itr->first << std::dec << "; Energy "
0206                                  << std::get<0>(itr->second) << "; Position (" << std::get<1>(itr->second).x() << ", "
0207                                  << std::get<1>(itr->second).y() << ", " << std::get<1>(itr->second).z() << ")";
0208     }
0209 #endif
0210   } else {
0211     edm::LogVerbatim("HGCalValid") << "No EE SimHit Found ";
0212   }
0213 
0214   //Accesing fh simhits
0215   const edm::Handle<std::vector<PCaloHit>>& fhSimHits = iEvent.getHandle(fhSimHitToken_);
0216   if (fhSimHits.isValid()) {
0217     analyzeHGCalSimHit(fhSimHits, 1, hefEnSim, fhHitRefs);
0218 #ifdef EDM_ML_DEBUG
0219     for (std::map<unsigned int, HGCHitTuple>::iterator itr = fhHitRefs.begin(); itr != fhHitRefs.end(); ++itr) {
0220       int idx = std::distance(fhHitRefs.begin(), itr);
0221       edm::LogInfo("HGCalValid") << "FHHit[" << idx << "] " << std::hex << itr->first << std::dec << "; Energy "
0222                                  << std::get<0>(itr->second) << "; Position (" << std::get<1>(itr->second).x() << ", "
0223                                  << std::get<1>(itr->second).y() << ", " << std::get<1>(itr->second).z() << ")";
0224     }
0225 #endif
0226   } else {
0227     edm::LogVerbatim("HGCalValid") << "No FH SimHit Found ";
0228   }
0229 
0230   //Accessing bh simhits
0231   const edm::Handle<std::vector<PCaloHit>>& bhSimHits = iEvent.getHandle(bhSimHitToken_);
0232   if (bhSimHits.isValid()) {
0233     analyzeHGCalSimHit(bhSimHits, 2, hebEnSim, bhHitRefs);
0234 #ifdef EDM_ML_DEBUG
0235     for (std::map<unsigned int, HGCHitTuple>::iterator itr = bhHitRefs.begin(); itr != bhHitRefs.end(); ++itr) {
0236       int idx = std::distance(bhHitRefs.begin(), itr);
0237       edm::LogInfo("HGCalValid") << "BHHit[" << idx << "] " << std::hex << itr->first << std::dec << "; Energy "
0238                                  << std::get<0>(itr->second) << "; Position (" << std::get<1>(itr->second).x() << ", "
0239                                  << std::get<1>(itr->second).y() << ", " << std::get<1>(itr->second).z() << ")";
0240     }
0241 #endif
0242   } else {
0243     edm::LogVerbatim("HGCalValid") << "No BH SimHit Found ";
0244   }
0245 
0246   //accessing EE Rechit information
0247   const edm::Handle<HGCeeRecHitCollection>& eeRecHit = iEvent.getHandle(eeRecHitToken_);
0248   if (eeRecHit.isValid()) {
0249     const HGCeeRecHitCollection* theHits = (eeRecHit.product());
0250     for (auto it = theHits->begin(); it != theHits->end(); ++it) {
0251       double energy = it->energy();
0252       heeEnRec->Fill(energy);
0253       std::map<unsigned int, HGCHitTuple>::const_iterator itr = eeHitRefs.find(it->id().rawId());
0254       if (itr != eeHitRefs.end()) {
0255         GlobalPoint xyz = hgcGeometry_[0]->getPosition(it->id());
0256         heeRecVsSimX->Fill(std::get<1>(itr->second).x(), xyz.x());
0257         heeRecVsSimY->Fill(std::get<1>(itr->second).y(), xyz.y());
0258         heeRecVsSimZ->Fill(std::get<1>(itr->second).z(), xyz.z());
0259         heedxVsX->Fill(std::get<1>(itr->second).x(), (xyz.x() - std::get<1>(itr->second).x()));
0260         heedyVsY->Fill(std::get<1>(itr->second).y(), (xyz.y() - std::get<1>(itr->second).y()));
0261         heedzVsZ->Fill(std::get<1>(itr->second).z(), (xyz.z() - std::get<1>(itr->second).z()));
0262         heeEnSimRec->Fill(std::get<0>(itr->second), energy);
0263 #ifdef EDM_ML_DEBUG
0264         edm::LogInfo("HGCalValid") << "EEHit: " << std::hex << it->id().rawId() << std::dec << " Sim ("
0265                                    << std::get<0>(itr->second) << ", " << std::get<1>(itr->second) << ") Rec ("
0266                                    << energy << ", " << xyz.x() << ", " << xyz.y() << ", " << xyz.z() << ")";
0267 #endif
0268       }
0269     }
0270   } else {
0271     edm::LogVerbatim("HGCalValid") << "No EE RecHit Found ";
0272   }
0273 
0274   //accessing FH Rechit information
0275   const edm::Handle<HGChefRecHitCollection>& fhRecHit = iEvent.getHandle(fhRecHitToken_);
0276   if (fhRecHit.isValid()) {
0277     const HGChefRecHitCollection* theHits = (fhRecHit.product());
0278     for (auto it = theHits->begin(); it != theHits->end(); ++it) {
0279       double energy = it->energy();
0280       hefEnRec->Fill(energy);
0281       std::map<unsigned int, HGCHitTuple>::const_iterator itr = fhHitRefs.find(it->id().rawId());
0282       if (itr != fhHitRefs.end()) {
0283         GlobalPoint xyz = hgcGeometry_[1]->getPosition(it->id());
0284 
0285         hefRecVsSimX->Fill(std::get<1>(itr->second).x(), xyz.x());
0286         hefRecVsSimY->Fill(std::get<1>(itr->second).y(), xyz.y());
0287         hefRecVsSimZ->Fill(std::get<1>(itr->second).z(), xyz.z());
0288         hefdxVsX->Fill(std::get<1>(itr->second).x(), (xyz.x() - std::get<1>(itr->second).x()));
0289         hefdyVsY->Fill(std::get<1>(itr->second).y(), (xyz.y() - std::get<1>(itr->second).y()));
0290         hefdzVsZ->Fill(std::get<1>(itr->second).z(), (xyz.z() - std::get<1>(itr->second).z()));
0291         hefEnSimRec->Fill(std::get<0>(itr->second), energy);
0292 #ifdef EDM_ML_DEBUG
0293         edm::LogInfo("HGCalValid") << "FHHit: " << std::hex << it->id().rawId() << std::dec << " Sim ("
0294                                    << std::get<0>(itr->second) << ", " << std::get<1>(itr->second) << ") Rec ("
0295                                    << energy << "," << xyz.x() << ", " << xyz.y() << ", " << xyz.z() << ")";
0296 #endif
0297       }
0298     }
0299   } else {
0300     edm::LogVerbatim("HGCalValid") << "No FH RecHit Found ";
0301   }
0302 
0303   //accessing BH Rechit information
0304   const edm::Handle<HGChebRecHitCollection>& bhRecHit = iEvent.getHandle(bhRecHitToken_);
0305   if (bhRecHit.isValid()) {
0306     const HGChebRecHitCollection* theHits = (bhRecHit.product());
0307     for (auto it = theHits->begin(); it != theHits->end(); ++it) {
0308       double energy = it->energy();
0309       hebEnRec->Fill(energy);
0310       std::map<unsigned int, HGCHitTuple>::const_iterator itr = bhHitRefs.find(it->id().rawId());
0311       GlobalPoint xyz = hgcGeometry_[2]->getPosition(it->id());
0312       if (itr != bhHitRefs.end()) {
0313         float ang3 = xyz.phi().value();  // returns the phi in radians
0314         float ang3s = std::get<1>(itr->second).phi().value();
0315         hebRecVsSimX->Fill(std::get<1>(itr->second).x(), xyz.x());
0316         hebRecVsSimY->Fill(std::get<1>(itr->second).y(), xyz.y());
0317         hebRecVsSimZ->Fill(std::get<1>(itr->second).z(), xyz.z());
0318         hebdEtaVsEta->Fill(std::get<1>(itr->second).eta(), (xyz.eta() - std::get<1>(itr->second).eta()));
0319         hebdPhiVsPhi->Fill(std::get<1>(itr->second).phi(), (ang3 - ang3s));
0320         hebdzVsZ->Fill(std::get<1>(itr->second).z(), (xyz.z() - std::get<1>(itr->second).z()));
0321         hebEnSimRec->Fill(std::get<0>(itr->second), energy);
0322 
0323 #ifdef EDM_ML_DEBUG
0324         edm::LogInfo("HGCalValid") << "BHHit: " << std::hex << it->id().rawId() << std::dec << " Sim ("
0325                                    << std::get<0>(itr->second) << ", " << std::get<1>(itr->second) << ") Rec ("
0326                                    << energy << ", " << xyz.x() << ", " << xyz.y() << ", " << xyz.z() << ")\n";
0327 #endif
0328       }
0329     }
0330   } else {
0331     edm::LogVerbatim("HGCalValid") << "No BH RecHit Found ";
0332   }
0333 }
0334 
0335 void HGCalHitValidation::analyzeHGCalSimHit(edm::Handle<std::vector<PCaloHit>> const& simHits,
0336                                             int idet,
0337                                             MonitorElement* hist,
0338                                             std::map<unsigned int, HGCHitTuple>& hitRefs) {
0339   for (std::vector<PCaloHit>::const_iterator simHit = simHits->begin(); simHit != simHits->end(); ++simHit) {
0340     DetId id(simHit->id());
0341     GlobalPoint p = hgcGeometry_[idet]->getPosition(id);
0342     float energy = simHit->energy();
0343 
0344     float energySum(energy);
0345     if (hitRefs.count(id.rawId()) != 0)
0346       energySum += std::get<0>(hitRefs[id.rawId()]);
0347     hitRefs[id.rawId()] = std::make_tuple(energySum, p);
0348     hist->Fill(energy);
0349   }
0350 }
0351 
0352 //define this as a plug-in
0353 #include "FWCore/Framework/interface/MakerMacros.h"
0354 DEFINE_FWK_MODULE(HGCalHitValidation);