Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-02-20 23:14:38

0001 #include "FWCore/Framework/interface/Frameworkfwd.h"
0002 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0003 
0004 #include "FWCore/Framework/interface/Event.h"
0005 #include "FWCore/Framework/interface/EventSetup.h"
0006 #include "FWCore/Framework/interface/MakerMacros.h"
0007 
0008 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0009 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0010 #include "FWCore/Utilities/interface/InputTag.h"
0011 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0012 #include "FWCore/Utilities/interface/Exception.h"
0013 
0014 #include "FWCore/ServiceRegistry/interface/Service.h"
0015 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0016 
0017 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0018 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0019 #include "DataFormats/Math/interface/Point3D.h"
0020 #include "DataFormats/Math/interface/Vector3D.h"
0021 
0022 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0023 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0024 #include "SimG4CMS/Calo/interface/CaloHitID.h"
0025 
0026 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0027 #include "Geometry/HcalCommonData/interface/HcalDDDRecConstants.h"
0028 #include "Geometry/HcalCommonData/interface/HcalHitRelabeller.h"
0029 #include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"
0030 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0031 #include "Geometry/Records/interface/HcalRecNumberingRecord.h"
0032 
0033 #include <TH2F.h>
0034 #include <cmath>
0035 #include <map>
0036 #include <memory>
0037 #include <string>
0038 #include <vector>
0039 
0040 class HcalSimHitAnalysis : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0041 public:
0042   HcalSimHitAnalysis(const edm::ParameterSet& ps);
0043   ~HcalSimHitAnalysis() override = default;
0044   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0045 
0046 protected:
0047   void beginJob() override;
0048   void analyze(edm::Event const&, edm::EventSetup const&) override;
0049   void beginRun(edm::Run const&, edm::EventSetup const&) override {}
0050   void endRun(edm::Run const&, edm::EventSetup const&) override {}
0051 
0052 private:
0053   static const int ndets_ = 4;
0054   const std::string g4Label_, hcalHits_;
0055   const bool verbose_, testNumber_;
0056   const edm::EDGetTokenT<edm::PCaloHitContainer> tok_calo_;
0057   const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> tok_geom_;
0058   const edm::ESGetToken<HcalDDDRecConstants, HcalRecNumberingRecord> tok_HRNDC_;
0059   TH2F *poszp_[ndets_], *poszn_[ndets_];
0060 };
0061 
0062 HcalSimHitAnalysis::HcalSimHitAnalysis(const edm::ParameterSet& ps)
0063     : g4Label_(ps.getUntrackedParameter<std::string>("moduleLabel", "g4SimHits")),
0064       hcalHits_(ps.getUntrackedParameter<std::string>("HitCollection", "HcalHits")),
0065       verbose_(ps.getUntrackedParameter<bool>("Verbose", false)),
0066       testNumber_(ps.getUntrackedParameter<bool>("TestNumber", true)),
0067       tok_calo_(consumes<edm::PCaloHitContainer>(edm::InputTag(g4Label_, hcalHits_))),
0068       tok_geom_(esConsumes<CaloGeometry, CaloGeometryRecord>()),
0069       tok_HRNDC_(esConsumes<HcalDDDRecConstants, HcalRecNumberingRecord>()) {
0070   usesResource(TFileService::kSharedResource);
0071 
0072   edm::LogVerbatim("HitStudy") << "Module Label: " << g4Label_ << "   Hits: " << hcalHits_ << " testNumber "
0073                                << testNumber_;
0074 }
0075 
0076 void HcalSimHitAnalysis::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0077   edm::ParameterSetDescription desc;
0078   desc.addUntracked<std::string>("ModuleLabel", "g4SimHits");
0079   desc.addUntracked<std::string>("HitCollection", "HcalHits");
0080   desc.addUntracked<bool>("Verbose", false);
0081   desc.addUntracked<bool>("TestNumber", true);
0082   descriptions.add("hcalSimHitAnalysis", desc);
0083 }
0084 
0085 void HcalSimHitAnalysis::beginJob() {
0086   edm::Service<TFileService> tfile;
0087   if (!tfile.isAvailable())
0088     throw cms::Exception("BadConfig") << "TFileService unavailable: "
0089                                       << "please add it to config file";
0090   char name[20], title[120];
0091   std::string dets[ndets_] = {"HB", "HE", "HO", "HF"};
0092   int nx[ndets_] = {100, 100, 350, 160};
0093   double xlo[ndets_] = {0, -300, 0, -160};
0094   double xhi[ndets_] = {500, 300, 3500, 160};
0095   int ny[ndets_] = {100, 100, 50, 160};
0096   double ylo[ndets_] = {170, -300, 375, -160};
0097   double yhi[ndets_] = {370, 300, 425, 160};
0098   std::string xttl[ndets_] = {"|z| (cm)", "x (cm)", "|z| (cm)", "x (cm)"};
0099   std::string yttl[ndets_] = {"#rho (cm)", "y (cm)", "#rho (cm)", "y (cm)"};
0100   for (int i = 0; i < ndets_; i++) {
0101     sprintf(name, "poszp%d", i);
0102     sprintf(title, "%s+", dets[i].c_str());
0103     poszp_[i] = tfile->make<TH2F>(name, title, nx[i], xlo[i], xhi[i], ny[i], ylo[i], yhi[i]);
0104     poszp_[i]->GetXaxis()->SetTitle(xttl[i].c_str());
0105     poszp_[i]->GetYaxis()->SetTitle(yttl[i].c_str());
0106     sprintf(title, "%s-", dets[i].c_str());
0107     poszp_[i]->GetYaxis()->SetTitleOffset(1.2);
0108     poszp_[i]->Sumw2();
0109     sprintf(name, "poszn%d", i);
0110     poszn_[i] = tfile->make<TH2F>(name, title, nx[i], xlo[i], xhi[i], ny[i], ylo[i], yhi[i]);
0111     poszn_[i]->GetXaxis()->SetTitle(xttl[i].c_str());
0112     poszn_[i]->GetYaxis()->SetTitle(yttl[i].c_str());
0113     poszn_[i]->GetYaxis()->SetTitleOffset(1.2);
0114     poszn_[i]->Sumw2();
0115   }
0116 }
0117 
0118 void HcalSimHitAnalysis::analyze(const edm::Event& e, const edm::EventSetup& iS) {
0119   if (verbose_)
0120     edm::LogVerbatim("HitStudy") << "Run = " << e.id().run() << " Event = " << e.id().event();
0121 
0122   // get hcalGeometry
0123   const CaloGeometry* geo = &iS.getData(tok_geom_);
0124   const HcalGeometry* hgeom = static_cast<const HcalGeometry*>(geo->getSubdetectorGeometry(DetId::Hcal, HcalBarrel));
0125   const auto& pHRNDC = iS.getData(tok_HRNDC_);
0126   const HcalDDDRecConstants* hcons = &pHRNDC;
0127 
0128   const edm::Handle<edm::PCaloHitContainer>& hitsCalo = e.getHandle(tok_calo_);
0129   bool getHits = (hitsCalo.isValid());
0130   uint32_t nhits = (getHits) ? hitsCalo->size() : 0;
0131   if (verbose_)
0132     edm::LogVerbatim("HitStudy") << "HcalSimHitAnalysis: Input flags Hits " << getHits << " with " << nhits << " hits";
0133   if (getHits) {
0134     std::vector<PCaloHit> hits;
0135     hits.insert(hits.end(), hitsCalo->begin(), hitsCalo->end());
0136     if (!hits.empty()) {
0137       std::map<HcalDetId, double> hitMap;
0138       for (auto hit : hits) {
0139         double edep = hit.energy();
0140         uint32_t id = hit.id();
0141         HcalDetId hid = (testNumber_) ? HcalHitRelabeller::relabel(id, hcons) : HcalDetId(id);
0142         auto it = hitMap.find(hid);
0143         if (it == hitMap.end()) {
0144           hitMap[hid] = edep;
0145         } else {
0146           (it->second) += edep;
0147         }
0148       }
0149 
0150       for (auto it : hitMap) {
0151         HcalDetId id(it.first);
0152         GlobalPoint gpos = hgeom->getPosition(id);
0153         HcalSubdetector subdet = (id).subdet();
0154         int indx =
0155             ((subdet == HcalBarrel)
0156                  ? 0
0157                  : ((subdet == HcalEndcap) ? 1 : ((subdet == HcalOuter) ? 2 : ((subdet == HcalForward) ? 3 : -1))));
0158         if (verbose_)
0159           edm::LogVerbatim("HitStudy") << "HcalSimHitAnalysis: " << id << ":" << it.second << " at " << gpos
0160                                        << " subdet " << subdet << ":" << indx;
0161         if (indx >= 0) {
0162           double x = ((indx == 0) || (indx == 2)) ? std::abs(gpos.z()) : gpos.x();
0163           double y = ((indx == 0) || (indx == 2)) ? (gpos.perp()) : gpos.y();
0164           if (id.zside() >= 0)
0165             poszp_[indx]->Fill(x, y);
0166           else
0167             poszn_[indx]->Fill(x, y);
0168         }
0169       }
0170     }
0171   }
0172 }
0173 
0174 //define this as a plug-in
0175 DEFINE_FWK_MODULE(HcalSimHitAnalysis);