HcalHits/plugins/HcalGeomCheck.cc

0001 // system include files
0002 #include <cmath>
0003 #include <iostream>
0004 #include <fstream>
0005 #include <vector>
0006 #include <map>
0007 #include <string>
0008
0009 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0010
0011 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0012
0013 #include "FWCore/Framework/interface/Frameworkfwd.h"
0014 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0015 #include "FWCore/Framework/interface/Event.h"
0016 #include "FWCore/Framework/interface/EventSetup.h"
0017 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0018 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0019 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0020 #include "FWCore/ServiceRegistry/interface/Service.h"
0021 #include "FWCore/Utilities/interface/InputTag.h"
0022
0023 #include "Geometry/HcalCommonData/interface/HcalHitRelabeller.h"
0024 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0025 #include "Geometry/HcalCommonData/interface/HcalDDDRecConstants.h"
0026 #include "Geometry/Records/interface/HcalRecNumberingRecord.h"
0027
0028 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0029 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0030
0031 #include <CLHEP/Geometry/Point3D.h>
0032 #include <CLHEP/Geometry/Transform3D.h>
0033 #include <CLHEP/Geometry/Vector3D.h>
0034 #include <CLHEP/Units/SystemOfUnits.h>
0035 #include <CLHEP/Units/PhysicalConstants.h>
0036
0037 #include "TH1D.h"
0038 #include "TH2D.h"
0039
0040 class HcalGeomCheck : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0041 public:
0042   struct hitsinfo {
0043     hitsinfo() { phi = eta = energy = time = 0.0; }
0044     double phi, eta, energy, time;
0045   };
0046
0047   explicit HcalGeomCheck(const edm::ParameterSet&);
0048   ~HcalGeomCheck() override = default;
0049   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0050
0051 protected:
0052   void beginJob() override;
0053   void beginRun(edm::Run const&, edm::EventSetup const&) override;
0054   void analyze(edm::Event const&, edm::EventSetup const&) override;
0055   void endRun(edm::Run const&, edm::EventSetup const&) override {}
0056
0057 private:
0058   void analyzeHits(int, const std::string&, const std::vector<PCaloHit>&);
0059
0060   // ----------member data ---------------------------
0061   const std::string caloHitSource_;
0062   const int ietaMin_, ietaMax_, depthMax_;
0063   const double rmin_, rmax_, zmin_, zmax_;
0064   const int nbinR_, nbinZ_, verbosity_;
0065   const edm::EDGetTokenT<edm::PCaloHitContainer> tok_hits_;
0066   const edm::ESGetToken<HcalDDDRecConstants, HcalRecNumberingRecord> tok_HRNDC_;
0067   const HcalDDDRecConstants* hcons_;
0068
0069   //histogram related stuff
0070   TH2D* h_RZ_;
0071   std::map<std::pair<int, int>, TH1D*> h_phi_;
0072   std::map<int, TH1D*> h_eta_;
0073   std::vector<TH1D*> h_E_, h_T_;
0074 };
0075
0076 HcalGeomCheck::HcalGeomCheck(const edm::ParameterSet& iConfig)
0077     : caloHitSource_(iConfig.getParameter<std::string>("caloHitSource")),
0078       ietaMin_(iConfig.getUntrackedParameter<int>("ietaMin", -41)),
0079       ietaMax_(iConfig.getUntrackedParameter<int>("ietaMax", 41)),
0080       depthMax_(iConfig.getUntrackedParameter<int>("depthMax", 7)),
0081       rmin_(iConfig.getUntrackedParameter<double>("rMin", 0.0)),
0082       rmax_(iConfig.getUntrackedParameter<double>("rMax", 5500.0)),
0083       zmin_(iConfig.getUntrackedParameter<double>("zMin", -12500.0)),
0084       zmax_(iConfig.getUntrackedParameter<double>("zMax", 12500.0)),
0085       nbinR_(iConfig.getUntrackedParameter<int>("nBinR", 550)),
0086       nbinZ_(iConfig.getUntrackedParameter<int>("nBinZ", 2500)),
0087       verbosity_(iConfig.getUntrackedParameter<int>("verbosity", 0)),
0088       tok_hits_(consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits", caloHitSource_))),
0089       tok_HRNDC_(esConsumes<HcalDDDRecConstants, HcalRecNumberingRecord, edm::Transition::BeginRun>()) {
0090   usesResource(TFileService::kSharedResource);
0091 }
0092
0093 void HcalGeomCheck::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0094   edm::ParameterSetDescription desc;
0095   desc.add<std::string>("caloHitSource", "HcalHits");
0096   desc.addUntracked<int>("ietaMin", -41);
0097   desc.addUntracked<int>("ietaMax", 41);
0098   desc.addUntracked<int>("depthMax", 7);
0099   desc.addUntracked<double>("rMin", 0.0);
0100   desc.addUntracked<double>("rMax", 5500.0);
0101   desc.addUntracked<double>("zMin", -12500.0);
0102   desc.addUntracked<double>("zMax", 12500.0);
0103   desc.addUntracked<int>("nBinR", 550);
0104   desc.addUntracked<int>("nBinZ", 250);
0105   desc.addUntracked<int>("verbosity", 0);
0106   descriptions.add("hcalGeomCheck", desc);
0107 }
0108
0109 void HcalGeomCheck::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0110   //Now the hits
0111   const edm::Handle<edm::PCaloHitContainer>& theCaloHitContainer = iEvent.getHandle(tok_hits_);
0112   if (theCaloHitContainer.isValid()) {
0113     if (verbosity_ > 0)
0114       edm::LogVerbatim("HcalValidation") << " PcalohitItr = " << theCaloHitContainer->size();
0115
0116     //Merge hits for the same DetID
0117     std::map<HcalDetId, hitsinfo> map_hits;
0118     for (auto const& hit : *(theCaloHitContainer.product())) {
0119       unsigned int id = hit.id();
0120       HcalDetId detId = HcalHitRelabeller::relabel(id, hcons_);
0121       double energy = hit.energy();
0122       double time = hit.time();
0123       int subdet = detId.subdet();
0124       int ieta = detId.ieta();
0125       int iphi = detId.iphi();
0126       int depth = detId.depth();
0127       std::pair<double, double> etaphi = hcons_->getEtaPhi(subdet, ieta, iphi);
0128       double rz = hcons_->getRZ(subdet, ieta, depth);
0129       if (verbosity_ > 2)
0130         edm::LogVerbatim("HcalValidation") << "i/p " << subdet << ":" << ieta << ":" << iphi << ":" << depth << " o/p "
0131                                            << etaphi.first << ":" << etaphi.second << ":" << rz;
0132       HepGeom::Point3D<float> gcoord = HepGeom::Point3D<float>(rz * cos(etaphi.second) / cosh(etaphi.first),
0133                                                                rz * sin(etaphi.second) / cosh(etaphi.first),
0134                                                                rz * tanh(etaphi.first));
0135       double tof = (gcoord.mag() * CLHEP::mm) / CLHEP::c_light;
0136       if (verbosity_ > 1)
0137         edm::LogVerbatim("HcalValidation")
0138             << "Detector " << subdet << " ieta = " << ieta << " iphi = " << iphi << " depth = " << depth
0139             << " positon = " << gcoord << " energy = " << energy << " time = " << time << ":" << tof;
0140       time -= tof;
0141       if (time < 0)
0142         time = 0;
0143       hitsinfo hinfo;
0144       if (map_hits.count(detId) != 0) {
0145         hinfo = map_hits[detId];
0146       } else {
0147         hinfo.phi = gcoord.getPhi();
0148         hinfo.eta = gcoord.getEta();
0149         hinfo.time = time;
0150       }
0151       hinfo.energy += energy;
0152       map_hits[detId] = hinfo;
0153
0154       h_RZ_->Fill(gcoord.z(), gcoord.rho());
0155     }
0156
0157     //Fill in histograms
0158     for (auto const& hit : map_hits) {
0159       hitsinfo hinfo = hit.second;
0160       int subdet = hit.first.subdet();
0161       if (subdet > 0 && subdet < static_cast<int>(h_E_.size())) {
0162         int depth = hit.first.depth();
0163         int ieta = hit.first.ieta();
0164         int iphi = hit.first.iphi();
0165         if (verbosity_ > 1)
0166           edm::LogVerbatim("HGCalValidation")
0167               << " ----------------------   eta = " << ieta << ":" << hinfo.eta << " phi = " << iphi << ":" << hinfo.phi
0168               << " depth = " << depth << " E = " << hinfo.energy << " T = " << hinfo.time;
0169         h_E_[subdet]->Fill(hinfo.energy);
0170         h_T_[subdet]->Fill(hinfo.time);
0171         auto itr1 = h_phi_.find(std::pair<int, int>(ieta, depth));
0172         if (itr1 != h_phi_.end())
0173           itr1->second->Fill(iphi);
0174         auto itr2 = h_eta_.find(depth);
0175         if (itr2 != h_eta_.end())
0176           itr2->second->Fill(ieta);
0177       }
0178     }
0179   } else if (verbosity_ > 0) {
0180     edm::LogVerbatim("HcalValidation") << "PCaloHitContainer does not "
0181                                        << "exist for HCAL";
0182   }
0183 }
0184
0185 // ------------ method called when starting to processes a run  ------------
0186 void HcalGeomCheck::beginRun(const edm::Run&, const edm::EventSetup& iSetup) {
0187   const auto& pHRNDC = iSetup.getData(tok_HRNDC_);
0188   hcons_ = &pHRNDC;
0189   if (verbosity_ > 0)
0190     edm::LogVerbatim("HcalValidation") << "Obtain HcalDDDRecConstants from Event Setup";
0191 }
0192
0193 void HcalGeomCheck::beginJob() {
0194   if (verbosity_ > 2)
0195     edm::LogVerbatim("HcalValidation") << "HcalGeomCheck:: Enter beginJob";
0196   edm::Service<TFileService> fs;
0197   h_RZ_ = fs->make<TH2D>("RZ", "R vs Z", nbinZ_, zmin_, zmax_, nbinR_, rmin_, rmax_);
0198   if (verbosity_ > 2)
0199     edm::LogVerbatim("HcalValidation") << "HcalGeomCheck: booked scatterplot RZ";
0200   char name[20], title[100];
0201   for (int depth = 1; depth < depthMax_; ++depth) {
0202     for (int ieta = ietaMin_; ieta <= ietaMax_; ++ieta) {
0203       sprintf(name, "phi%d%d", ieta, depth);
0204       sprintf(title, "i#phi (i#eta = %d, depth = %d)", ieta, depth);
0205       h_phi_[std::pair<int, int>(ieta, depth)] = fs->make<TH1D>(name, title, 400, -20, 380);
0206       if (verbosity_ > 2)
0207         edm::LogVerbatim("HcalValidation") << "HcalGeomCheck: books " << title;
0208     }
0209     sprintf(name, "eta%d", depth);
0210     sprintf(title, "i#eta (depth = %d)", depth);
0211     h_eta_[depth] = fs->make<TH1D>(name, title, 100, -50, 50);
0212     if (verbosity_ > 2)
0213       edm::LogVerbatim("HcalValidation") << "HcalGeomCheck: books " << title;
0214   }
0215
0216   std::vector<std::string> dets = {"HCAL", "HB", "HE", "HF", "HO"};
0217   for (unsigned int ih = 0; ih < dets.size(); ++ih) {
0218     sprintf(name, "E_%s", dets[ih].c_str());
0219     sprintf(title, "Energy deposit in %s (MeV)", dets[ih].c_str());
0220     h_E_.emplace_back(fs->make<TH1D>(name, title, 1000, 0.0, 1.0));
0221     if (verbosity_ > 2)
0222       edm::LogVerbatim("HcalValidation") << "HcalGeomCheck: books " << title;
0223     sprintf(name, "T_%s", dets[ih].c_str());
0224     sprintf(title, "Time of hit in %s (ns)", dets[ih].c_str());
0225     h_T_.emplace_back(fs->make<TH1D>(name, title, 1000, 0.0, 200.0));
0226     if (verbosity_ > 2)
0227       edm::LogVerbatim("HcalValidation") << "HcalGeomCheck: books " << title;
0228   }
0229 }
0230
0231 #include "FWCore/Framework/interface/MakerMacros.h"
0232 //define this as a plug-in
0233 DEFINE_FWK_MODULE(HcalGeomCheck);