Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-03-11 02:41:11

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/GlobalSystemOfUnits.h"
0035 #include "CLHEP/Units/GlobalPhysicalConstants.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     unsigned int nused(0);
0119     for (auto const& hit : *(theCaloHitContainer.product())) {
0120       unsigned int id = hit.id();
0121       HcalDetId detId = HcalHitRelabeller::relabel(id, hcons_);
0122       double energy = hit.energy();
0123       double time = hit.time();
0124       int subdet = detId.subdet();
0125       int ieta = detId.ieta();
0126       int iphi = detId.iphi();
0127       int depth = detId.depth();
0128       std::pair<double, double> etaphi = hcons_->getEtaPhi(subdet, ieta, iphi);
0129       double rz = hcons_->getRZ(subdet, ieta, depth);
0130       if (verbosity_ > 2)
0131         edm::LogVerbatim("HcalValidation") << "i/p " << subdet << ":" << ieta << ":" << iphi << ":" << depth << " o/p "
0132                                            << etaphi.first << ":" << etaphi.second << ":" << rz;
0133       HepGeom::Point3D<float> gcoord = HepGeom::Point3D<float>(rz * cos(etaphi.second) / cosh(etaphi.first),
0134                                                                rz * sin(etaphi.second) / cosh(etaphi.first),
0135                                                                rz * tanh(etaphi.first));
0136       nused++;
0137       double tof = (gcoord.mag() * CLHEP::mm) / CLHEP::c_light;
0138       if (verbosity_ > 1)
0139         edm::LogVerbatim("HcalValidation")
0140             << "Detector " << subdet << " ieta = " << ieta << " iphi = " << iphi << " depth = " << depth
0141             << " positon = " << gcoord << " energy = " << energy << " time = " << time << ":" << tof;
0142       time -= tof;
0143       if (time < 0)
0144         time = 0;
0145       hitsinfo hinfo;
0146       if (map_hits.count(detId) != 0) {
0147         hinfo = map_hits[detId];
0148       } else {
0149         hinfo.phi = gcoord.getPhi();
0150         hinfo.eta = gcoord.getEta();
0151         hinfo.time = time;
0152       }
0153       hinfo.energy += energy;
0154       map_hits[detId] = hinfo;
0155 
0156       h_RZ_->Fill(gcoord.z(), gcoord.rho());
0157     }
0158 
0159     //Fill in histograms
0160     for (auto const& hit : map_hits) {
0161       hitsinfo hinfo = hit.second;
0162       int subdet = hit.first.subdet();
0163       if (subdet > 0 && subdet < static_cast<int>(h_E_.size())) {
0164         int depth = hit.first.depth();
0165         int ieta = hit.first.ieta();
0166         int iphi = hit.first.iphi();
0167         if (verbosity_ > 1)
0168           edm::LogVerbatim("HGCalValidation")
0169               << " ----------------------   eta = " << ieta << ":" << hinfo.eta << " phi = " << iphi << ":" << hinfo.phi
0170               << " depth = " << depth << " E = " << hinfo.energy << " T = " << hinfo.time;
0171         h_E_[subdet]->Fill(hinfo.energy);
0172         h_T_[subdet]->Fill(hinfo.time);
0173         auto itr1 = h_phi_.find(std::pair<int, int>(ieta, depth));
0174         if (itr1 != h_phi_.end())
0175           itr1->second->Fill(iphi);
0176         auto itr2 = h_eta_.find(depth);
0177         if (itr2 != h_eta_.end())
0178           itr2->second->Fill(ieta);
0179       }
0180     }
0181   } else if (verbosity_ > 0) {
0182     edm::LogVerbatim("HcalValidation") << "PCaloHitContainer does not "
0183                                        << "exist for HCAL";
0184   }
0185 }
0186 
0187 // ------------ method called when starting to processes a run  ------------
0188 void HcalGeomCheck::beginRun(const edm::Run&, const edm::EventSetup& iSetup) {
0189   const auto& pHRNDC = iSetup.getData(tok_HRNDC_);
0190   hcons_ = &pHRNDC;
0191   if (verbosity_ > 0)
0192     edm::LogVerbatim("HcalValidation") << "Obtain HcalDDDRecConstants from Event Setup";
0193 }
0194 
0195 void HcalGeomCheck::beginJob() {
0196   if (verbosity_ > 2)
0197     edm::LogVerbatim("HcalValidation") << "HcalGeomCheck:: Enter beginJob";
0198   edm::Service<TFileService> fs;
0199   h_RZ_ = fs->make<TH2D>("RZ", "R vs Z", nbinZ_, zmin_, zmax_, nbinR_, rmin_, rmax_);
0200   if (verbosity_ > 2)
0201     edm::LogVerbatim("HcalValidation") << "HcalGeomCheck: booked scatterplot RZ";
0202   char name[20], title[100];
0203   for (int depth = 1; depth < depthMax_; ++depth) {
0204     for (int ieta = ietaMin_; ieta <= ietaMax_; ++ieta) {
0205       sprintf(name, "phi%d%d", ieta, depth);
0206       sprintf(title, "i#phi (i#eta = %d, depth = %d)", ieta, depth);
0207       h_phi_[std::pair<int, int>(ieta, depth)] = fs->make<TH1D>(name, title, 400, -20, 380);
0208       if (verbosity_ > 2)
0209         edm::LogVerbatim("HcalValidation") << "HcalGeomCheck: books " << title;
0210     }
0211     sprintf(name, "eta%d", depth);
0212     sprintf(title, "i#eta (depth = %d)", depth);
0213     h_eta_[depth] = fs->make<TH1D>(name, title, 100, -50, 50);
0214     if (verbosity_ > 2)
0215       edm::LogVerbatim("HcalValidation") << "HcalGeomCheck: books " << title;
0216   }
0217 
0218   std::vector<std::string> dets = {"HCAL", "HB", "HE", "HF", "HO"};
0219   for (unsigned int ih = 0; ih < dets.size(); ++ih) {
0220     sprintf(name, "E_%s", dets[ih].c_str());
0221     sprintf(title, "Energy deposit in %s (MeV)", dets[ih].c_str());
0222     h_E_.emplace_back(fs->make<TH1D>(name, title, 1000, 0.0, 1.0));
0223     if (verbosity_ > 2)
0224       edm::LogVerbatim("HcalValidation") << "HcalGeomCheck: books " << title;
0225     sprintf(name, "T_%s", dets[ih].c_str());
0226     sprintf(title, "Time of hit in %s (ns)", dets[ih].c_str());
0227     h_T_.emplace_back(fs->make<TH1D>(name, title, 1000, 0.0, 200.0));
0228     if (verbosity_ > 2)
0229       edm::LogVerbatim("HcalValidation") << "HcalGeomCheck: books " << title;
0230   }
0231 }
0232 
0233 #include "FWCore/Framework/interface/MakerMacros.h"
0234 //define this as a plug-in
0235 DEFINE_FWK_MODULE(HcalGeomCheck);