File indexing completed on 2022-12-03 01:10:33
0001
0002 #include <iostream>
0003 #include <fstream>
0004 #include <vector>
0005 #include <map>
0006 #include <string>
0007
0008
0009 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0010 #include "DataFormats/HGCDigi/interface/HGCDigiCollections.h"
0011
0012 #include "FWCore/Framework/interface/Frameworkfwd.h"
0013 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0014 #include "FWCore/Framework/interface/Event.h"
0015 #include "FWCore/Framework/interface/EventSetup.h"
0016 #include "FWCore/Framework/interface/MakerMacros.h"
0017 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0018 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0019 #include "FWCore/ServiceRegistry/interface/Service.h"
0020
0021 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0022 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0023 #include "DataFormats/ForwardDetId/interface/HGCScintillatorDetId.h"
0024
0025
0026 #include "TROOT.h"
0027 #include "TSystem.h"
0028 #include "TFile.h"
0029 #include "TH1.h"
0030 #include "TH2.h"
0031
0032 class HGCalBHValidation : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0033 public:
0034 explicit HGCalBHValidation(const edm::ParameterSet& ps);
0035 ~HGCalBHValidation() override = default;
0036
0037 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0038
0039 protected:
0040 void beginJob() override {}
0041 void beginRun(edm::Run const&, edm::EventSetup const&) override;
0042 void endRun(edm::Run const&, edm::EventSetup const&) override {}
0043 void analyze(edm::Event const&, edm::EventSetup const&) override;
0044
0045 private:
0046 edm::Service<TFileService> fs_;
0047 const std::string g4Label_, hits_;
0048 const edm::InputTag digis_;
0049 const int iSample_;
0050 const double threshold_;
0051 const edm::EDGetTokenT<edm::PCaloHitContainer> tok_hits_;
0052 const edm::EDGetTokenT<HGCalDigiCollection> tok_digi_;
0053 const int ringMax_;
0054
0055 TH1D *hsimE1_, *hsimE2_, *hsimTm_;
0056 TH1D *hsimLn_, *hdigEn_, *hdigLn_;
0057 TH2D *hsimOc_, *hsi2Oc_, *hsi3Oc_;
0058 TH2D *hdigOc_, *hdi2Oc_, *hdi3Oc_;
0059 };
0060
0061 HGCalBHValidation::HGCalBHValidation(const edm::ParameterSet& ps)
0062 : g4Label_(ps.getParameter<std::string>("ModuleLabel")),
0063 hits_((ps.getParameter<std::string>("HitCollection"))),
0064 digis_(ps.getParameter<edm::InputTag>("DigiCollection")),
0065 iSample_(ps.getParameter<int>("Sample")),
0066 threshold_(ps.getParameter<double>("Threshold")),
0067 tok_hits_(consumes<edm::PCaloHitContainer>(edm::InputTag(g4Label_, hits_))),
0068 tok_digi_(consumes<HGCalDigiCollection>(digis_)),
0069 ringMax_(50) {
0070 usesResource(TFileService::kSharedResource);
0071
0072 edm::LogVerbatim("HGCalValidation") << "HGCalBHValidation::Input for SimHit:" << edm::InputTag(g4Label_, hits_)
0073 << " Digits:" << digis_ << " Sample: " << iSample_ << " Threshold "
0074 << threshold_;
0075 }
0076
0077 void HGCalBHValidation::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0078 edm::ParameterSetDescription desc;
0079 desc.add<std::string>("ModuleLabel", "g4SimHits");
0080 desc.add<std::string>("HitCollection", "HGCHitsHEback");
0081 desc.add<edm::InputTag>("DigiCollection", edm::InputTag("simHGCalUnsuppressedDigis", "HEback"));
0082 desc.add<int>("Sample", 5);
0083 desc.add<double>("Threshold", 15.0);
0084 descriptions.add("hgcalBHAnalysis", desc);
0085 }
0086
0087 void HGCalBHValidation::beginRun(edm::Run const&, edm::EventSetup const& es) {
0088 edm::LogVerbatim("HGCalValidation") << "HGCalBHValidation::Maximum Number of"
0089 << " ring sectors:" << ringMax_ << "\nHitsValidationHcal::Booking the Histograms";
0090
0091
0092 hsimE1_ = fs_->make<TH1D>("SimHitEn1", "Sim Hit Energy", 1000, 0.0, 1.0);
0093 hsimE2_ = fs_->make<TH1D>("SimHitEn2", "Sim Hit Energy", 1000, 0.0, 1.0);
0094 hsimTm_ = fs_->make<TH1D>("SimHitTime", "Sim Hit Time", 1000, 0.0, 500.0);
0095 hsimLn_ = fs_->make<TH1D>("SimHitLong", "Sim Hit Long. Profile", 50, 0.0, 25.0);
0096 hsimOc_ = fs_->make<TH2D>("SimHitOccup", "Sim Hit Occupnacy", 2 * ringMax_ + 1, -ringMax_, ringMax_ + 1, 360, 0, 360);
0097 hsi2Oc_ = fs_->make<TH2D>("SimHitOccu2", "Sim Hit Occupnacy", 2 * ringMax_ + 1, -ringMax_, ringMax_ + 1, 360, 0, 360);
0098 hsi3Oc_ = fs_->make<TH2D>("SimHitOccu3", "Sim Hit Occupnacy", 2 * ringMax_ + 1, -ringMax_, ringMax_ + 1, 50, 0, 25);
0099
0100 hdigEn_ = fs_->make<TH1D>("DigiEnergy", "Digi ADC Sample", 1000, 0.0, 1000.0);
0101 hdigLn_ = fs_->make<TH1D>("DigiLong", "Digi Long. Profile", 50, 0.0, 25.0);
0102 hdigOc_ = fs_->make<TH2D>("DigiOccup", "Digi Occupnacy", 2 * ringMax_ + 1, -ringMax_, ringMax_ + 1, 360, 0, 360);
0103 hdi2Oc_ = fs_->make<TH2D>("DigiOccu2", "Digi Occupnacy", 2 * ringMax_ + 1, -ringMax_, ringMax_ + 1, 360, 0, 360);
0104 hdi3Oc_ = fs_->make<TH2D>("DigiOccu3", "Digi Occupnacy", 2 * ringMax_ + 1, -ringMax_, ringMax_ + 1, 50, 0, 25);
0105 }
0106
0107 void HGCalBHValidation::analyze(const edm::Event& e, const edm::EventSetup&) {
0108 edm::LogVerbatim("HGCalValidation") << "HGCalBHValidation:Run = " << e.id().run() << " Event = " << e.id().event();
0109
0110
0111 const edm::Handle<edm::PCaloHitContainer>& hitsHE = e.getHandle(tok_hits_);
0112 edm::LogVerbatim("HGCalValidation") << "HGCalBHValidation.: PCaloHitContainer"
0113 << " obtained with flag " << hitsHE.isValid();
0114 if (hitsHE.isValid()) {
0115 edm::LogVerbatim("HGCalValidation") << "HGCalBHValidation: PCaloHit buffer " << hitsHE->size();
0116 unsigned i(0);
0117 std::map<unsigned int, double> map_try;
0118 for (edm::PCaloHitContainer::const_iterator it = hitsHE->begin(); it != hitsHE->end(); ++it) {
0119 double energy = it->energy();
0120 double time = it->time();
0121 unsigned int id = it->id();
0122 if (DetId(id).det() == DetId::HGCalHSc) {
0123 double ring = 0.01 + HGCScintillatorDetId(id).ring();
0124 double phi = HGCScintillatorDetId(id).iphi() - 0.01;
0125 int lay = HGCScintillatorDetId(id).layer();
0126 double ring1 = ring * HGCScintillatorDetId(id).zside();
0127 hsi2Oc_->Fill(ring1, phi, energy);
0128 hsimE1_->Fill(energy);
0129 hsimTm_->Fill(time, energy);
0130 hsimOc_->Fill(ring1, phi, energy);
0131 hsimLn_->Fill(lay, energy);
0132 hsi3Oc_->Fill(ring1, lay, energy);
0133 double ensum(0);
0134 if (map_try.count(id) != 0)
0135 ensum = map_try[id];
0136 ensum += energy;
0137 map_try[id] = ensum;
0138 edm::LogVerbatim("HGCalValidation") << "HGCalBHHit[" << ++i << "] ID " << std::hex << " " << id << std::dec
0139 << " " << HGCScintillatorDetId(id) << " E " << energy << " time " << time;
0140 }
0141 }
0142 for (std::map<unsigned int, double>::iterator itr = map_try.begin(); itr != map_try.end(); ++itr) {
0143 hsimE2_->Fill((*itr).second);
0144 }
0145 }
0146
0147
0148 unsigned int kount(0);
0149 const edm::Handle<HGCalDigiCollection>& hecoll = e.getHandle(tok_digi_);
0150 edm::LogVerbatim("HGCalValidation") << "HGCalBHValidation.: "
0151 << "HGCalDigiCollection obtained with"
0152 << " flag " << hecoll.isValid();
0153 if (hecoll.isValid()) {
0154 edm::LogVerbatim("HGCalValidation") << "HGCalBHValidation: HGCalDigi "
0155 << "buffer " << hecoll->size();
0156 for (HGCalDigiCollection::const_iterator it = hecoll->begin(); it != hecoll->end(); ++it) {
0157 const HGCalDataFrame& df(*it);
0158 double energy = df[iSample_].data();
0159 if (DetId(df.id()).det() == DetId::HGCalHSc) {
0160 HGCScintillatorDetId cell(df.id());
0161 int depth = cell.layer();
0162 if (energy > threshold_) {
0163 double ring = 0.01 + cell.ring();
0164 double phi = cell.iphi() - 0.01;
0165 double ring1 = cell.zside() * ring;
0166 hdi2Oc_->Fill(ring1, phi);
0167 hdigEn_->Fill(energy);
0168 hdigOc_->Fill(ring1, phi);
0169 hdigLn_->Fill(depth);
0170 hdi3Oc_->Fill(ring1, depth);
0171 edm::LogVerbatim("HGCalValidation")
0172 << "HGCalBHDigit[" << ++kount << "] ID " << cell << " E " << energy << ":" << (energy > threshold_);
0173 }
0174 }
0175 }
0176 }
0177 }
0178
0179 DEFINE_FWK_MODULE(HGCalBHValidation);
