HGCalValidation/test/HGCHitValidation.cc

0001 /// -*- C++ -*-
0002 //
0003 // Package:    HGCHitValidation
0004 // Class:      HGCHitValidation
0005 //
0006 /**\class HGCHitValidation HGCHitValidation.cc Validation/HGCalValidation/test/HGCHitValidation.cc
0007
0008  Description: [one line class summary]
0009
0010  Implementation:
0011     [Notes on implementation]
0012 */
0013 //
0014 // Original Author:  "Maksat Haytmyradov"
0015 //         Created:  Fri March 19 13:32:26 CDT 2016
0016 // $Id$
0017 //
0018 //
0019
0020 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0021
0022 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0023 #include "Geometry/HGCalCommonData/interface/HGCalGeometryMode.h"
0024 #include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
0025 #include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h"
0026
0027 #include "DataFormats/Common/interface/Handle.h"
0028 #include "DataFormats/HGCRecHit/interface/HGCRecHit.h"
0029 #include "DataFormats/HGCRecHit/interface/HGCRecHitCollections.h"
0030 #include "DataFormats/ForwardDetId/interface/HGCalDetId.h"
0031 #include "DataFormats/ForwardDetId/interface/HGCSiliconDetId.h"
0032 #include "DataFormats/ForwardDetId/interface/HGCScintillatorDetId.h"
0033 #include "DataFormats/ForwardDetId/interface/ForwardSubdetector.h"
0034
0035 #include "FWCore/Framework/interface/Frameworkfwd.h"
0036 #include "FWCore/Framework/interface/Event.h"
0037 #include "FWCore/Framework/interface/MakerMacros.h"
0038 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0039 #include "FWCore/Framework/interface/ESHandle.h"
0040 #include "FWCore/ServiceRegistry/interface/Service.h"
0041 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0042 #include "FWCore/Utilities/interface/Exception.h"
0043 #include "FWCore/Utilities/interface/EDGetToken.h"
0044 #include "FWCore/Utilities/interface/transform.h"
0045
0046 #include "SimG4CMS/Calo/interface/HGCNumberingScheme.h"
0047 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0048 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0049 #include "SimDataFormats/CaloTest/interface/HGCalTestNumbering.h"
0050
0051 #include <TH1.h>
0052 #include <TH2.h>
0053 #include <TTree.h>
0054 #include <TVector3.h>
0055 #include <TSystem.h>
0056 #include <TFile.h>
0057
0058 #include <cmath>
0059 #include <memory>
0060 #include <iostream>
0061 #include <string>
0062 #include <vector>
0063
0064 class HGCHitValidation : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0065 public:
0066   explicit HGCHitValidation(const edm::ParameterSet &);
0067   ~HGCHitValidation() override = default;
0068   static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0069
0070 private:
0071   typedef std::tuple<float, float, float, float> HGCHitTuple;
0072
0073   void beginJob() override;
0074   void endJob() override {}
0075   void beginRun(edm::Run const &, edm::EventSetup const &) override;
0076   void analyze(edm::Event const &, edm::EventSetup const &) override;
0077   void endRun(edm::Run const &, edm::EventSetup const &) override {}
0078   virtual void beginLuminosityBlock(edm::LuminosityBlock const &, edm::EventSetup const &) {}
0079   virtual void endLuminosityBlock(edm::LuminosityBlock const &, edm::EventSetup const &) {}
0080   void analyzeHGCalSimHit(edm::Handle<std::vector<PCaloHit>> const &simHits,
0081                           int idet,
0082                           TH1F *,
0083                           std::map<unsigned int, HGCHitTuple> &);
0084   template <class T1>
0085   void analyzeHGCalRecHit(T1 const &theHits, std::map<unsigned int, HGCHitTuple> const &hitRefs);
0086
0087 private:
0088   //HGC Geometry
0089   const std::vector<std::string> geometrySource_;
0090   const std::vector<int> ietaExcludeBH_;
0091   const bool makeTree_;
0092   const std::vector<edm::ESGetToken<HGCalDDDConstants, IdealGeometryRecord>> tok_hgcal_;
0093   const std::vector<edm::ESGetToken<HGCalGeometry, IdealGeometryRecord>> tok_hgcalg_;
0094   std::vector<const HGCalDDDConstants *> hgcCons_;
0095   std::vector<const HGCalGeometry *> hgcGeometry_;
0096
0097   const edm::InputTag eeSimHitSource, fhSimHitSource, bhSimHitSource;
0098   const edm::EDGetTokenT<std::vector<PCaloHit>> eeSimHitToken_;
0099   const edm::EDGetTokenT<std::vector<PCaloHit>> fhSimHitToken_;
0100   const edm::EDGetTokenT<std::vector<PCaloHit>> bhSimHitToken_;
0101   const edm::InputTag eeRecHitSource, fhRecHitSource, bhRecHitSource;
0102   const edm::EDGetTokenT<HGCeeRecHitCollection> eeRecHitToken_;
0103   const edm::EDGetTokenT<HGChefRecHitCollection> fhRecHitToken_;
0104   const edm::EDGetTokenT<HGChebRecHitCollection> bhRecHitToken_;
0105
0106   TTree *hgcHits_;
0107   std::vector<float> *heeRecX_, *heeRecY_, *heeRecZ_, *heeRecEnergy_;
0108   std::vector<float> *hefRecX_, *hefRecY_, *hefRecZ_, *hefRecEnergy_;
0109   std::vector<float> *hebRecX_, *hebRecY_, *hebRecZ_, *hebRecEnergy_;
0110   std::vector<float> *heeSimX_, *heeSimY_, *heeSimZ_, *heeSimEnergy_;
0111   std::vector<float> *hefSimX_, *hefSimY_, *hefSimZ_, *hefSimEnergy_;
0112   std::vector<float> *hebSimX_, *hebSimY_, *hebSimZ_, *hebSimEnergy_;
0113   std::vector<float> *hebSimEta_, *hebRecEta_, *hebSimPhi_, *hebRecPhi_;
0114   std::vector<unsigned int> *heeDetID_, *hefDetID_, *hebDetID_;
0115
0116   TH2F *heedzVsZ_, *heedyVsY_, *heedxVsX_;
0117   TH2F *hefdzVsZ_, *hefdyVsY_, *hefdxVsX_;
0118   TH2F *hebdzVsZ_, *hebdPhiVsPhi_, *hebdEtaVsEta_;
0119   TH2F *heeRecVsSimZ_, *heeRecVsSimY_, *heeRecVsSimX_;
0120   TH2F *hefRecVsSimZ_, *hefRecVsSimY_, *hefRecVsSimX_;
0121   TH2F *hebRecVsSimZ_, *hebRecVsSimY_, *hebRecVsSimX_;
0122   TH2F *heeEnSimRec_, *hefEnSimRec_, *hebEnSimRec_;
0123   TH1F *hebEnRec_, *hebEnSim_, *hefEnRec_;
0124   TH1F *hefEnSim_, *heeEnRec_, *heeEnSim_;
0125 };
0126
0127 HGCHitValidation::HGCHitValidation(const edm::ParameterSet &cfg)
0128     : geometrySource_(cfg.getUntrackedParameter<std::vector<std::string>>("geometrySource")),
0129       ietaExcludeBH_(cfg.getParameter<std::vector<int>>("ietaExcludeBH")),
0130       makeTree_(cfg.getUntrackedParameter<bool>("makeTree", true)),
0131       tok_hgcal_{
0132           edm::vector_transform(geometrySource_,
0133                                 [this](const std::string &name) {
0134                                   return esConsumes<HGCalDDDConstants, IdealGeometryRecord, edm::Transition::BeginRun>(
0135                                       edm::ESInputTag{"", name});
0136                                 })},
0137       tok_hgcalg_{edm::vector_transform(
0138           geometrySource_,
0139           [this](const std::string &name) {
0140             return esConsumes<HGCalGeometry, IdealGeometryRecord, edm::Transition::BeginRun>(edm::ESInputTag{"", name});
0141           })},
0142       eeSimHitSource(cfg.getParameter<edm::InputTag>("eeSimHitSource")),
0143       fhSimHitSource(cfg.getParameter<edm::InputTag>("fhSimHitSource")),
0144       bhSimHitSource(cfg.getParameter<edm::InputTag>("bhSimHitSource")),
0145       eeSimHitToken_(consumes<std::vector<PCaloHit>>(eeSimHitSource)),
0146       fhSimHitToken_(consumes<std::vector<PCaloHit>>(fhSimHitSource)),
0147       bhSimHitToken_(consumes<std::vector<PCaloHit>>(bhSimHitSource)),
0148       eeRecHitSource(cfg.getParameter<edm::InputTag>("eeRecHitSource")),
0149       fhRecHitSource(cfg.getParameter<edm::InputTag>("fhRecHitSource")),
0150       bhRecHitSource(cfg.getParameter<edm::InputTag>("bhRecHitSource")),
0151       eeRecHitToken_(consumes<HGCeeRecHitCollection>(eeRecHitSource)),
0152       fhRecHitToken_(consumes<HGChefRecHitCollection>(fhRecHitSource)),
0153       bhRecHitToken_(consumes<HGChebRecHitCollection>(bhRecHitSource)) {
0154   usesResource(TFileService::kSharedResource);
0155
0156   hgcHits_ = nullptr;
0157   heeRecX_ = heeRecY_ = heeRecZ_ = heeRecEnergy_ = nullptr;
0158   hefRecX_ = hefRecY_ = hefRecZ_ = hefRecEnergy_ = nullptr;
0159   hebRecX_ = hebRecY_ = hebRecZ_ = hebRecEnergy_ = nullptr;
0160   heeSimX_ = heeSimY_ = heeSimZ_ = heeSimEnergy_ = nullptr;
0161   hefSimX_ = hefSimY_ = hefSimZ_ = hefSimEnergy_ = nullptr;
0162   hebSimX_ = hebSimY_ = hebSimZ_ = hebSimEnergy_ = nullptr;
0163   hebSimEta_ = hebRecEta_ = hebSimPhi_ = hebRecPhi_ = nullptr;
0164   heeDetID_ = hefDetID_ = hebDetID_ = nullptr;
0165   heedzVsZ_ = heedyVsY_ = heedxVsX_ = nullptr;
0166   hefdzVsZ_ = hefdyVsY_ = hefdxVsX_ = nullptr;
0167   hebdzVsZ_ = hebdPhiVsPhi_ = hebdEtaVsEta_ = nullptr;
0168   heeRecVsSimZ_ = heeRecVsSimY_ = heeRecVsSimX_ = nullptr;
0169   hefRecVsSimZ_ = hefRecVsSimY_ = hefRecVsSimX_ = nullptr;
0170   hebRecVsSimZ_ = hebRecVsSimY_ = hebRecVsSimX_ = nullptr;
0171   heeEnSimRec_ = hefEnSimRec_ = hebEnSimRec_ = nullptr;
0172   hebEnRec_ = hebEnSim_ = hefEnRec_ = nullptr;
0173   hefEnSim_ = heeEnRec_ = heeEnSim_ = nullptr;
0174
0175   edm::LogVerbatim("HGCalValid") << "MakeTree Flag set to " << makeTree_ << " and use " << geometrySource_.size()
0176                                  << " Geometry sources";
0177   for (auto const &s : geometrySource_)
0178     edm::LogVerbatim("HGCalValid") << "  " << s;
0179   edm::LogVerbatim("HGCalValid") << "SimHit labels: " << eeSimHitSource << "  " << fhSimHitSource << "  "
0180                                  << bhSimHitSource;
0181   edm::LogVerbatim("HGCalValid") << "RecHit labels: " << eeRecHitSource << "  " << fhRecHitSource << "  "
0182                                  << bhRecHitSource;
0183   edm::LogVerbatim("HGCalValid") << "Exclude the following " << ietaExcludeBH_.size() << " ieta values from BH plots";
0184   for (unsigned int k = 0; k < ietaExcludeBH_.size(); ++k)
0185     edm::LogVerbatim("HGCalValid") << " [" << k << "] " << ietaExcludeBH_[k];
0186 }
0187
0188 void HGCHitValidation::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0189   std::vector<std::string> names = {"HGCalEESensitive", "HGCalHESiliconSensitive", "HGCalHEScintillatorSensitive"};
0190   std::vector<int> etas;
0191   edm::ParameterSetDescription desc;
0192   desc.addUntracked<bool>("makeTree", true);
0193   desc.addUntracked<std::vector<std::string>>("geometrySource", names);
0194   desc.add<edm::InputTag>("eeSimHitSource", edm::InputTag("g4SimHits", "HGCHitsEE"));
0195   desc.add<edm::InputTag>("fhSimHitSource", edm::InputTag("g4SimHits", "HGCHitsHEfront"));
0196   desc.add<edm::InputTag>("bhSimHitSource", edm::InputTag("g4SimHits", "HGCHitsHEback"));
0197   desc.add<edm::InputTag>("eeRecHitSource", edm::InputTag("HGCalRecHit", "HGCEERecHits"));
0198   desc.add<edm::InputTag>("fhRecHitSource", edm::InputTag("HGCalRecHit", "HGCHEFRecHits"));
0199   desc.add<edm::InputTag>("bhRecHitSource", edm::InputTag("HGCalRecHit", "HGCHEBRecHits"));
0200   desc.add<std::vector<int>>("ietaExcludeBH", etas);
0201   descriptions.add("hgcHitAnalysis", desc);
0202 }
0203
0204 void HGCHitValidation::beginJob() {
0205   //initiating fileservice
0206   edm::Service<TFileService> fs;
0207   if (makeTree_) {
0208     hgcHits_ = fs->make<TTree>("hgcHits", "Hit Collection");
0209     hgcHits_->Branch("heeRecX", &heeRecX_);
0210     hgcHits_->Branch("heeRecY", &heeRecY_);
0211     hgcHits_->Branch("heeRecZ", &heeRecZ_);
0212     hgcHits_->Branch("heeRecEnergy", &heeRecEnergy_);
0213     hgcHits_->Branch("hefRecX", &hefRecX_);
0214     hgcHits_->Branch("hefRecY", &hefRecY_);
0215     hgcHits_->Branch("hefRecZ", &hefRecZ_);
0216     hgcHits_->Branch("hefRecEnergy", &hefRecEnergy_);
0217     hgcHits_->Branch("hebRecX", &hebRecX_);
0218     hgcHits_->Branch("hebRecY", &hebRecY_);
0219     hgcHits_->Branch("hebRecZ", &hebRecZ_);
0220     hgcHits_->Branch("hebRecEta", &hebRecEta_);
0221     hgcHits_->Branch("hebRecPhi", &hebRecPhi_);
0222     hgcHits_->Branch("hebRecEnergy", &hebRecEnergy_);
0223
0224     hgcHits_->Branch("heeSimX", &heeSimX_);
0225     hgcHits_->Branch("heeSimY", &heeSimY_);
0226     hgcHits_->Branch("heeSimZ", &heeSimZ_);
0227     hgcHits_->Branch("heeSimEnergy", &heeSimEnergy_);
0228     hgcHits_->Branch("hefSimX", &hefSimX_);
0229     hgcHits_->Branch("hefSimY", &hefSimY_);
0230     hgcHits_->Branch("hefSimZ", &hefSimZ_);
0231     hgcHits_->Branch("hefSimEnergy", &hefSimEnergy_);
0232     hgcHits_->Branch("hebSimX", &hebSimX_);
0233     hgcHits_->Branch("hebSimY", &hebSimY_);
0234     hgcHits_->Branch("hebSimZ", &hebSimZ_);
0235     hgcHits_->Branch("hebSimEta", &hebSimEta_);
0236     hgcHits_->Branch("hebSimPhi", &hebSimPhi_);
0237     hgcHits_->Branch("hebSimEnergy", &hebSimEnergy_);
0238
0239     hgcHits_->Branch("heeDetID", &heeDetID_);
0240     hgcHits_->Branch("hefDetID", &hefDetID_);
0241     hgcHits_->Branch("hebDetID", &hebDetID_);
0242   } else {
0243     heedzVsZ_ = fs->make<TH2F>("heedzVsZ", "", 7200, -360, 360, 100, -0.1, 0.1);
0244     heedyVsY_ = fs->make<TH2F>("heedyVsY", "", 400, -200, 200, 100, -0.02, 0.02);
0245     heedxVsX_ = fs->make<TH2F>("heedxVsX", "", 400, -200, 200, 100, -0.02, 0.02);
0246     heeRecVsSimZ_ = fs->make<TH2F>("heeRecVsSimZ", "", 7200, -360, 360, 7200, -360, 360);
0247     heeRecVsSimY_ = fs->make<TH2F>("heeRecVsSimY", "", 400, -200, 200, 400, -200, 200);
0248     heeRecVsSimX_ = fs->make<TH2F>("heeRecVsSimX", "", 400, -200, 200, 400, -200, 200);
0249     hefdzVsZ_ = fs->make<TH2F>("hefdzVsZ", "", 8200, -410, 410, 100, -0.1, 0.1);
0250     hefdyVsY_ = fs->make<TH2F>("hefdyVsY", "", 400, -200, 200, 100, -0.02, 0.02);
0251     hefdxVsX_ = fs->make<TH2F>("hefdxVsX", "", 400, -200, 200, 100, -0.02, 0.02);
0252     hefRecVsSimZ_ = fs->make<TH2F>("hefRecVsSimZ", "", 8200, -410, 410, 8200, -410, 410);
0253     hefRecVsSimY_ = fs->make<TH2F>("hefRecVsSimY", "", 400, -200, 200, 400, -200, 200);
0254     hefRecVsSimX_ = fs->make<TH2F>("hefRecVsSimX", "", 400, -200, 200, 400, -200, 200);
0255     hebdzVsZ_ = fs->make<TH2F>("hebdzVsZ", "", 1080, -540, 540, 100, -1.0, 1.0);
0256     hebdPhiVsPhi_ = fs->make<TH2F>("hebdPhiVsPhi", "", M_PI * 100, -0.5, M_PI + 0.5, 200, -0.2, 0.2);
0257     hebdEtaVsEta_ = fs->make<TH2F>("hebdEtaVsEta", "", 1000, -5, 5, 200, -0.1, 0.1);
0258     hebRecVsSimZ_ = fs->make<TH2F>("hebRecVsSimZ", "", 1080, -540, 540, 1080, -540, 540);
0259     hebRecVsSimY_ = fs->make<TH2F>("hebRecVsSimY", "", 400, -200, 200, 400, -200, 200);
0260     hebRecVsSimX_ = fs->make<TH2F>("hebRecVsSimX", "", 400, -200, 200, 400, -200, 200);
0261     heeEnRec_ = fs->make<TH1F>("heeEnRec", "", 1000, 0, 10);
0262     heeEnSim_ = fs->make<TH1F>("heeEnSim", "", 1000, 0, 0.01);
0263     heeEnSimRec_ = fs->make<TH2F>("heeEnSimRec", "", 200, 0, 0.002, 200, 0, 0.2);
0264     hefEnRec_ = fs->make<TH1F>("hefEnRec", "", 1000, 0, 10);
0265     hefEnSim_ = fs->make<TH1F>("hefEnSim", "", 1000, 0, 0.01);
0266     hefEnSimRec_ = fs->make<TH2F>("hefEnSimRec", "", 200, 0, 0.001, 200, 0, 0.5);
0267     hebEnRec_ = fs->make<TH1F>("hebEnRec", "", 1000, 0, 15);
0268     hebEnSim_ = fs->make<TH1F>("hebEnSim", "", 1000, 0, 0.01);
0269     hebEnSimRec_ = fs->make<TH2F>("hebEnSimRec", "", 200, 0, 0.02, 200, 0, 4);
0270   }
0271 }
0272
0273 void HGCHitValidation::beginRun(edm::Run const &iRun, edm::EventSetup const &iSetup) {
0274   //initiating hgc Geometry
0275   for (size_t i = 0; i < geometrySource_.size(); i++) {
0276     edm::LogVerbatim("HGCalValid") << "Tries to initialize HGCalGeometry and HGCalDDDConstants for " << i;
0277     edm::ESHandle<HGCalDDDConstants> hgcCons = iSetup.getHandle(tok_hgcal_[i]);
0278     if (hgcCons.isValid()) {
0279       hgcCons_.push_back(hgcCons.product());
0280     } else {
0281       edm::LogWarning("HGCalValid") << "Cannot initiate HGCalDDDConstants for " << geometrySource_[i] << std::endl;
0282     }
0283     edm::ESHandle<HGCalGeometry> hgcGeom = iSetup.getHandle(tok_hgcalg_[i]);
0284     if (hgcGeom.isValid()) {
0285       hgcGeometry_.push_back(hgcGeom.product());
0286     } else {
0287       edm::LogWarning("HGCalValid") << "Cannot initiate HGCalGeometry for " << geometrySource_[i] << std::endl;
0288     }
0289   }
0290 }
0291
0292 void HGCHitValidation::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0293   std::map<unsigned int, HGCHitTuple> eeHitRefs, fhHitRefs, bhHitRefs;
0294
0295   //Accesing ee simhits
0296   const edm::Handle<std::vector<PCaloHit>> &eeSimHits = iEvent.getHandle(eeSimHitToken_);
0297
0298   if (eeSimHits.isValid()) {
0299     analyzeHGCalSimHit(eeSimHits, 0, heeEnSim_, eeHitRefs);
0300     for (std::map<unsigned int, HGCHitTuple>::iterator itr = eeHitRefs.begin(); itr != eeHitRefs.end(); ++itr) {
0301       int idx = std::distance(eeHitRefs.begin(), itr);
0302       edm::LogVerbatim("HGCalValid") << "EEHit[" << idx << "] " << std::hex << itr->first << std::dec << "; Energy "
0303                                      << std::get<0>(itr->second) << "; Position"
0304                                      << " (" << std::get<1>(itr->second) << ", " << std::get<2>(itr->second) << ", "
0305                                      << std::get<3>(itr->second) << ")";
0306     }
0307   } else {
0308     edm::LogWarning("HGCalValid") << "No EE SimHit Found " << std::endl;
0309   }
0310
0311   //Accesing fh simhits
0312   const edm::Handle<std::vector<PCaloHit>> &fhSimHits = iEvent.getHandle(fhSimHitToken_);
0313   if (fhSimHits.isValid()) {
0314     analyzeHGCalSimHit(fhSimHits, 1, hefEnSim_, fhHitRefs);
0315     for (std::map<unsigned int, HGCHitTuple>::iterator itr = fhHitRefs.begin(); itr != fhHitRefs.end(); ++itr) {
0316       int idx = std::distance(fhHitRefs.begin(), itr);
0317       edm::LogVerbatim("HGCalValid") << "FHHit[" << idx << "] " << std::hex << itr->first << std::dec << "; Energy "
0318                                      << std::get<0>(itr->second) << "; Position"
0319                                      << " (" << std::get<1>(itr->second) << ", " << std::get<2>(itr->second) << ", "
0320                                      << std::get<3>(itr->second) << ")";
0321     }
0322   } else {
0323     edm::LogWarning("HGCalValid") << "No FH SimHit Found " << std::endl;
0324   }
0325
0326   //Accessing bh simhits
0327   const edm::Handle<std::vector<PCaloHit>> &bhSimHits = iEvent.getHandle(bhSimHitToken_);
0328   if (bhSimHits.isValid()) {
0329     analyzeHGCalSimHit(bhSimHits, 2, hebEnSim_, bhHitRefs);
0330     for (std::map<unsigned int, HGCHitTuple>::iterator itr = bhHitRefs.begin(); itr != bhHitRefs.end(); ++itr) {
0331       int idx = std::distance(bhHitRefs.begin(), itr);
0332       edm::LogVerbatim("HGCalValid") << "BHHit[" << idx << "] " << std::hex << itr->first << std::dec << "; Energy "
0333                                      << std::get<0>(itr->second) << "; Position (" << std::get<1>(itr->second) << ", "
0334                                      << std::get<2>(itr->second) << ", " << std::get<3>(itr->second) << ")";
0335     }
0336   } else {
0337     edm::LogWarning("HGCalValid") << "No BH SimHit Found " << std::endl;
0338   }
0339
0340   //accessing EE Rechit information
0341   const edm::Handle<HGCeeRecHitCollection> &eeRecHit = iEvent.getHandle(eeRecHitToken_);
0342   if (eeRecHit.isValid()) {
0343     const HGCeeRecHitCollection *theHits = (eeRecHit.product());
0344     for (auto it = theHits->begin(); it != theHits->end(); ++it) {
0345       double energy = it->energy();
0346       if (!makeTree_)
0347         heeEnRec_->Fill(energy);
0348       std::map<unsigned int, HGCHitTuple>::const_iterator itr = eeHitRefs.find(it->id().rawId());
0349       if (itr != eeHitRefs.end()) {
0350         GlobalPoint xyz = hgcGeometry_[0]->getPosition(it->id());
0351         if (makeTree_) {
0352           heeRecX_->push_back(xyz.x());
0353           heeRecY_->push_back(xyz.y());
0354           heeRecZ_->push_back(xyz.z());
0355           heeRecEnergy_->push_back(energy);
0356           heeSimX_->push_back(std::get<1>(itr->second));
0357           heeSimY_->push_back(std::get<2>(itr->second));
0358           heeSimZ_->push_back(std::get<3>(itr->second));
0359           heeSimEnergy_->push_back(std::get<0>(itr->second));
0360           heeDetID_->push_back(itr->first);
0361         } else {
0362           heeRecVsSimX_->Fill(std::get<1>(itr->second), xyz.x());
0363           heeRecVsSimY_->Fill(std::get<2>(itr->second), xyz.y());
0364           heeRecVsSimZ_->Fill(std::get<3>(itr->second), xyz.z());
0365           heedxVsX_->Fill(std::get<1>(itr->second), (xyz.x() - std::get<1>(itr->second)));
0366           heedyVsY_->Fill(std::get<2>(itr->second), (xyz.y() - std::get<2>(itr->second)));
0367           heedzVsZ_->Fill(std::get<3>(itr->second), (xyz.z() - std::get<3>(itr->second)));
0368           heeEnSimRec_->Fill(std::get<0>(itr->second), energy);
0369         }
0370         edm::LogVerbatim("HGCalValid") << "EEHit: " << std::hex << it->id().rawId() << std::dec << " Sim ("
0371                                        << std::get<0>(itr->second) << ", " << std::get<1>(itr->second) << ", "
0372                                        << std::get<2>(itr->second) << ", " << std::get<3>(itr->second) << ") Rec ("
0373                                        << energy << ", " << xyz.x() << ", " << xyz.y() << ", " << xyz.z();
0374       }
0375     }
0376   } else {
0377     edm::LogWarning("HGCalValid") << "No EE RecHit Found " << std::endl;
0378   }
0379
0380   //accessing FH Rechit information
0381   const edm::Handle<HGChefRecHitCollection> &fhRecHit = iEvent.getHandle(fhRecHitToken_);
0382   if (fhRecHit.isValid()) {
0383     const HGChefRecHitCollection *theHits = (fhRecHit.product());
0384     for (auto it = theHits->begin(); it != theHits->end(); ++it) {
0385       double energy = it->energy();
0386       if (!makeTree_)
0387         hefEnRec_->Fill(energy);
0388       std::map<unsigned int, HGCHitTuple>::const_iterator itr = fhHitRefs.find(it->id().rawId());
0389       if (itr != fhHitRefs.end()) {
0390         GlobalPoint xyz = hgcGeometry_[1]->getPosition(it->id());
0391         if (makeTree_) {
0392           hefRecX_->push_back(xyz.x());
0393           hefRecY_->push_back(xyz.y());
0394           hefRecZ_->push_back(xyz.z());
0395           hefRecEnergy_->push_back(energy);
0396           hefSimX_->push_back(std::get<1>(itr->second));
0397           hefSimY_->push_back(std::get<2>(itr->second));
0398           hefSimZ_->push_back(std::get<3>(itr->second));
0399           hefSimEnergy_->push_back(std::get<0>(itr->second));
0400           hefDetID_->push_back(itr->first);
0401         } else {
0402           hefRecVsSimX_->Fill(std::get<1>(itr->second), xyz.x());
0403           hefRecVsSimY_->Fill(std::get<2>(itr->second), xyz.y());
0404           hefRecVsSimZ_->Fill(std::get<3>(itr->second), xyz.z());
0405           hefdxVsX_->Fill(std::get<1>(itr->second), (xyz.x() - std::get<1>(itr->second)));
0406           hefdyVsY_->Fill(std::get<2>(itr->second), (xyz.y() - std::get<2>(itr->second)));
0407           hefdzVsZ_->Fill(std::get<3>(itr->second), (xyz.z() - std::get<3>(itr->second)));
0408           hefEnSimRec_->Fill(std::get<0>(itr->second), energy);
0409         }
0410         edm::LogVerbatim("HGCalValid") << "FHHit: " << std::hex << it->id().rawId() << std::dec << " Sim ("
0411                                        << std::get<0>(itr->second) << ", " << std::get<1>(itr->second) << ", "
0412                                        << std::get<2>(itr->second) << ", " << std::get<3>(itr->second) << ") Rec ("
0413                                        << energy << "," << xyz.x() << ", " << xyz.y() << ", " << xyz.z();
0414       }
0415     }
0416   } else {
0417     edm::LogWarning("HGCalValid") << "No FH RecHit Found " << std::endl;
0418   }
0419
0420   //accessing BH Rechit information
0421   const edm::Handle<HGChebRecHitCollection> &bhRecHit = iEvent.getHandle(bhRecHitToken_);
0422   if (bhRecHit.isValid()) {
0423     const HGChebRecHitCollection *theHits = (bhRecHit.product());
0424     analyzeHGCalRecHit(theHits, bhHitRefs);
0425   } else {
0426     edm::LogWarning("HGCalValid") << "No BH RecHit Found " << std::endl;
0427   }
0428
0429   if (makeTree_) {
0430     hgcHits_->Fill();
0431
0432     heeRecX_->clear();
0433     heeRecY_->clear();
0434     heeRecZ_->clear();
0435     hefRecX_->clear();
0436     hefRecY_->clear();
0437     hefRecZ_->clear();
0438     hebRecX_->clear();
0439     hebRecY_->clear();
0440     hebRecZ_->clear();
0441     heeRecEnergy_->clear();
0442     hefRecEnergy_->clear();
0443     hebRecEnergy_->clear();
0444     heeSimX_->clear();
0445     heeSimY_->clear();
0446     heeSimZ_->clear();
0447     hefSimX_->clear();
0448     hefSimY_->clear();
0449     hefSimZ_->clear();
0450     hebSimX_->clear();
0451     hebSimY_->clear();
0452     hebSimZ_->clear();
0453     heeSimEnergy_->clear();
0454     hefSimEnergy_->clear();
0455     hebSimEnergy_->clear();
0456     hebSimEta_->clear();
0457     hebRecEta_->clear();
0458     hebSimPhi_->clear();
0459     hebRecPhi_->clear();
0460     heeDetID_->clear();
0461     hefDetID_->clear();
0462     hebDetID_->clear();
0463   }
0464 }
0465
0466 void HGCHitValidation::analyzeHGCalSimHit(edm::Handle<std::vector<PCaloHit>> const &simHits,
0467                                           int idet,
0468                                           TH1F *hist,
0469                                           std::map<unsigned int, HGCHitTuple> &hitRefs) {
0470   const HGCalTopology &hTopo = hgcGeometry_[idet]->topology();
0471   for (auto const &simHit : *simHits) {
0472     unsigned int id = simHit.id();
0473     std::pair<float, float> xy;
0474     bool ok(true);
0475     int subdet(0), zside, layer, wafer, celltype, cell, wafer2(0), cell2(0);
0476     if (hgcCons_[idet]->waferHexagon8()) {
0477       HGCSiliconDetId detId = HGCSiliconDetId(id);
0478       subdet = (int)(detId.det());
0479       cell = detId.cellU();
0480       cell2 = detId.cellV();
0481       wafer = detId.waferU();
0482       wafer2 = detId.waferV();
0483       celltype = detId.type();
0484       layer = detId.layer();
0485       zside = detId.zside();
0486       xy = hgcCons_[idet]->locateCell(zside, layer, wafer, wafer2, cell, cell2, false, true, false, false, false);
0487     } else if (hgcCons_[idet]->tileTrapezoid()) {
0488       HGCScintillatorDetId detId = HGCScintillatorDetId(id);
0489       subdet = (int)(detId.det());
0490       cell = detId.ietaAbs();
0491       wafer = detId.iphi();
0492       celltype = detId.type();
0493       layer = detId.layer();
0494       zside = detId.zside();
0495       xy = hgcCons_[idet]->locateCellTrap(zside, layer, wafer, cell, false, false);
0496     } else {
0497       HGCalTestNumbering::unpackHexagonIndex(simHit.id(), subdet, zside, layer, wafer, celltype, cell);
0498       xy = hgcCons_[idet]->locateCell(cell, layer, wafer, false);
0499
0500       //skip this hit if after ganging it is not valid
0501       std::pair<int, int> recoLayerCell = hgcCons_[idet]->simToReco(cell, layer, wafer, hTopo.detectorType());
0502       ok = !(recoLayerCell.second < 0 || recoLayerCell.first < 0);
0503       id = HGCalDetId((ForwardSubdetector)(subdet), zside, layer, celltype, wafer, cell).rawId();
0504     }
0505
0506     edm::LogVerbatim("HGCalValid") << "SimHit: " << std::hex << id << std::dec << " (" << subdet << ":" << zside << ":"
0507                                    << layer << ":" << celltype << ":" << wafer << ":" << wafer2 << ":" << cell << ":"
0508                                    << cell2 << ") Flag " << ok;
0509
0510     if (ok) {
0511       float zp = hgcCons_[idet]->waferZ(layer, false);
0512       if (zside < 0)
0513         zp = -zp;
0514       float xp = (zside < 0) ? -xy.first / 10 : xy.first / 10;
0515       float yp = xy.second / 10.0;
0516       float energy = simHit.energy();
0517
0518       float energySum(energy);
0519       if (hitRefs.count(id) != 0)
0520         energySum += std::get<0>(hitRefs[id]);
0521       hitRefs[id] = std::make_tuple(energySum, xp, yp, zp);
0522       if (hist != nullptr)
0523         hist->Fill(energy);
0524       edm::LogVerbatim("HGCalValid") << "Position (" << xp << ", " << yp << ", " << zp << ") "
0525                                      << " Energy " << simHit.energy() << ":" << energySum;
0526     }
0527   }
0528 }
0529
0530 template <class T1>
0531 void HGCHitValidation::analyzeHGCalRecHit(T1 const &theHits, std::map<unsigned int, HGCHitTuple> const &hitRefs) {
0532   for (auto it = theHits->begin(); it != theHits->end(); ++it) {
0533     DetId id = it->id();
0534     double energy = it->energy();
0535     if (!makeTree_)
0536       hebEnRec_->Fill(energy);
0537     GlobalPoint xyz = hgcGeometry_[2]->getPosition(id);
0538
0539     std::map<unsigned int, HGCHitTuple>::const_iterator itr = hitRefs.find(id.rawId());
0540     if (itr != hitRefs.end()) {
0541       float ang3 = xyz.phi().value();  // returns the phi in radians
0542       double fac = sinh(std::get<1>(itr->second));
0543       double pT = std::get<3>(itr->second) / fac;
0544       double xp = pT * cos(std::get<2>(itr->second));
0545       double yp = pT * sin(std::get<2>(itr->second));
0546       if (makeTree_) {
0547         hebRecX_->push_back(xyz.x());
0548         hebRecY_->push_back(xyz.y());
0549         hebRecZ_->push_back(xyz.z());
0550         hebRecEnergy_->push_back(energy);
0551         hebSimX_->push_back(xp);
0552         hebSimY_->push_back(yp);
0553         hebSimZ_->push_back(std::get<3>(itr->second));
0554         hebSimEnergy_->push_back(std::get<0>(itr->second));
0555         hebSimEta_->push_back(std::get<1>(itr->second));
0556         hebRecEta_->push_back(xyz.eta());
0557         hebSimPhi_->push_back(std::get<2>(itr->second));
0558         hebRecPhi_->push_back(ang3);
0559         hebDetID_->push_back(itr->first);
0560       } else {
0561         hebRecVsSimX_->Fill(xp, xyz.x());
0562         hebRecVsSimY_->Fill(yp, xyz.y());
0563         hebRecVsSimZ_->Fill(std::get<3>(itr->second), xyz.z());
0564         hebdEtaVsEta_->Fill(std::get<1>(itr->second), (xyz.eta() - std::get<1>(itr->second)));
0565         hebdPhiVsPhi_->Fill(std::get<2>(itr->second), (ang3 - std::get<2>(itr->second)));
0566         hebdzVsZ_->Fill(std::get<3>(itr->second), (xyz.z() - std::get<3>(itr->second)));
0567         hebEnSimRec_->Fill(std::get<0>(itr->second), energy);
0568       }
0569       edm::LogVerbatim("HGCalValid") << "BHHit: " << std::hex << id.rawId() << std::dec << " Sim ("
0570                                      << std::get<0>(itr->second) << ", " << std::get<1>(itr->second) << ", "
0571                                      << std::get<2>(itr->second) << ", " << std::get<3>(itr->second) << ") Rec ("
0572                                      << energy << ", " << xyz.eta() << ", " << ang3 << ", " << xyz.z() << ")";
0573     }
0574   }
0575 }
0576
0577 //define this as a plug-in
0578 DEFINE_FWK_MODULE(HGCHitValidation);