File indexing completed on 2021-05-05 03:16:18
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 {}
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 template <class T>
0045 void analyzeDigi(const T&, double const&, bool const&, int const&, unsigned int&);
0046
0047 private:
0048 edm::Service<TFileService> fs_;
0049 const std::string g4Label_, hits_;
0050 const edm::InputTag digis_;
0051 const int iSample_;
0052 const double threshold_;
0053 const edm::EDGetTokenT<edm::PCaloHitContainer> tok_hits_;
0054 const edm::EDGetToken tok_digi_;
0055 const int etaMax_;
0056
0057 TH1D *hsimE1_, *hsimE2_, *hsimTm_;
0058 TH1D *hsimLn_, *hdigEn_, *hdigLn_;
0059 TH2D *hsimOc_, *hsi2Oc_, *hsi3Oc_;
0060 TH2D *hdigOc_, *hdi2Oc_, *hdi3Oc_;
0061 };
0062
0063 HGCalBHValidation::HGCalBHValidation(const edm::ParameterSet& ps)
0064 : g4Label_(ps.getParameter<std::string>("ModuleLabel")),
0065 hits_((ps.getParameter<std::string>("HitCollection"))),
0066 digis_(ps.getParameter<edm::InputTag>("DigiCollection")),
0067 iSample_(ps.getParameter<int>("Sample")),
0068 threshold_(ps.getParameter<double>("Threshold")),
0069 tok_hits_(consumes<edm::PCaloHitContainer>(edm::InputTag(g4Label_, hits_))),
0070 tok_digi_(consumes<HGCalDigiCollection>(digis_)),
0071 etaMax_(100) {
0072 usesResource(TFileService::kSharedResource);
0073
0074 edm::LogVerbatim("HGCalValidation") << "HGCalBHValidation::Input for SimHit:" << edm::InputTag(g4Label_, hits_)
0075 << " Digits:" << digis_ << " Sample: " << iSample_ << " Threshold "
0076 << threshold_;
0077 }
0078
0079 void HGCalBHValidation::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0080 edm::ParameterSetDescription desc;
0081 desc.add<std::string>("ModuleLabel", "g4SimHits");
0082 desc.add<std::string>("HitCollection", "HGCHitsHEback");
0083 desc.add<edm::InputTag>("DigiCollection", edm::InputTag("simHGCalUnsuppressedDigis", "HEback"));
0084 desc.add<int>("Sample", 5);
0085 desc.add<double>("Threshold", 15.0);
0086 descriptions.add("hgcalBHAnalysis", desc);
0087 }
0088
0089 void HGCalBHValidation::beginRun(edm::Run const&, edm::EventSetup const& es) {
0090 edm::LogVerbatim("HGCalValidation") << "HGCalBHValidation::Maximum Number of"
0091 << " eta sectors:" << etaMax_ << "\nHitsValidationHcal::Booking the Histograms";
0092
0093
0094 hsimE1_ = fs_->make<TH1D>("SimHitEn1", "Sim Hit Energy", 1000, 0.0, 1.0);
0095 hsimE2_ = fs_->make<TH1D>("SimHitEn2", "Sim Hit Energy", 1000, 0.0, 1.0);
0096 hsimTm_ = fs_->make<TH1D>("SimHitTime", "Sim Hit Time", 1000, 0.0, 500.0);
0097 hsimLn_ = fs_->make<TH1D>("SimHitLong", "Sim Hit Long. Profile", 50, 0.0, 25.0);
0098 hsimOc_ = fs_->make<TH2D>("SimHitOccup", "Sim Hit Occupnacy", 2 * etaMax_ + 1, -etaMax_, etaMax_ + 1, 360, 0, 360);
0099 hsi2Oc_ = fs_->make<TH2D>("SimHitOccu2", "Sim Hit Occupnacy", 2 * etaMax_ + 1, -etaMax_, etaMax_ + 1, 360, 0, 360);
0100 hsi3Oc_ = fs_->make<TH2D>("SimHitOccu3", "Sim Hit Occupnacy", 2 * etaMax_ + 1, -etaMax_, etaMax_ + 1, 50, 0, 25);
0101
0102 hdigEn_ = fs_->make<TH1D>("DigiEnergy", "Digi ADC Sample", 1000, 0.0, 1000.0);
0103 hdigLn_ = fs_->make<TH1D>("DigiLong", "Digi Long. Profile", 50, 0.0, 25.0);
0104 hdigOc_ = fs_->make<TH2D>("DigiOccup", "Digi Occupnacy", 2 * etaMax_ + 1, -etaMax_, etaMax_ + 1, 360, 0, 360);
0105 hdi2Oc_ = fs_->make<TH2D>("DigiOccu2", "Digi Occupnacy", 2 * etaMax_ + 1, -etaMax_, etaMax_ + 1, 360, 0, 360);
0106 hdi3Oc_ = fs_->make<TH2D>("DigiOccu3", "Digi Occupnacy", 2 * etaMax_ + 1, -etaMax_, etaMax_ + 1, 50, 0, 25);
0107 }
0108
0109 void HGCalBHValidation::analyze(const edm::Event& e, const edm::EventSetup&) {
0110 edm::LogVerbatim("HGCalValidation") << "HGCalBHValidation:Run = " << e.id().run() << " Event = " << e.id().event();
0111
0112
0113 edm::Handle<edm::PCaloHitContainer> hitsHE;
0114 e.getByToken(tok_hits_, hitsHE);
0115 edm::LogVerbatim("HGCalValidation") << "HGCalBHValidation.: PCaloHitContainer"
0116 << " obtained with flag " << hitsHE.isValid();
0117 if (hitsHE.isValid()) {
0118 edm::LogVerbatim("HGCalValidation") << "HGCalBHValidation: PCaloHit buffer " << hitsHE->size();
0119 unsigned i(0);
0120 std::map<unsigned int, double> map_try;
0121 for (edm::PCaloHitContainer::const_iterator it = hitsHE->begin(); it != hitsHE->end(); ++it) {
0122 double energy = it->energy();
0123 double time = it->time();
0124 unsigned int id = it->id();
0125 int eta(0), phi(0), lay(0);
0126 bool bh = (DetId(id).det() == DetId::HGCalHSc);
0127 if (bh) {
0128 eta = HGCScintillatorDetId(id).ieta();
0129 phi = HGCScintillatorDetId(id).iphi();
0130 lay = HGCScintillatorDetId(id).layer();
0131 }
0132 double eta1 = (eta >= 0) ? (eta + 0.1) : (eta - 0.1);
0133 if (bh) {
0134 hsi2Oc_->Fill(eta1, (phi - 0.1), energy);
0135 hsimE1_->Fill(energy);
0136 hsimTm_->Fill(time, energy);
0137 hsimOc_->Fill(eta1, (phi - 0.1), energy);
0138 hsimLn_->Fill(lay, energy);
0139 hsi3Oc_->Fill(eta1, lay, energy);
0140 double ensum(0);
0141 if (map_try.count(id) != 0)
0142 ensum = map_try[id];
0143 ensum += energy;
0144 map_try[id] = ensum;
0145 ++i;
0146 edm::LogVerbatim("HGCalValidation") << "HGCalBHHit[" << i << "] ID " << std::hex << " " << id << std::dec << " "
0147 << HGCScintillatorDetId(id) << " E " << energy << " time " << time;
0148 }
0149 }
0150 for (std::map<unsigned int, double>::iterator itr = map_try.begin(); itr != map_try.end(); ++itr) {
0151 hsimE2_->Fill((*itr).second);
0152 }
0153 }
0154
0155
0156 unsigned int kount(0);
0157 edm::Handle<HGCalDigiCollection> hecoll;
0158 e.getByToken(tok_digi_, hecoll);
0159 edm::LogVerbatim("HGCalValidation") << "HGCalBHValidation.: "
0160 << "HGCalDigiCollection obtained with"
0161 << " flag " << hecoll.isValid();
0162 if (hecoll.isValid()) {
0163 edm::LogVerbatim("HGCalValidation") << "HGCalBHValidation: HGCalDigi "
0164 << "buffer " << hecoll->size();
0165 for (HGCalDigiCollection::const_iterator it = hecoll->begin(); it != hecoll->end(); ++it) {
0166 HGCalDataFrame df(*it);
0167 double energy = df[iSample_].data();
0168 bool bh = (DetId(df.id()).det() == DetId::HGCalHSc);
0169 if (bh) {
0170 HGCScintillatorDetId cell(df.id());
0171 int depth = cell.layer();
0172 analyzeDigi(cell, energy, bh, depth, kount);
0173 }
0174 }
0175 }
0176 }
0177
0178 template <class T>
0179 void HGCalBHValidation::analyzeDigi(
0180 const T& cell, double const& energy, bool const& bh, int const& depth, unsigned int& kount) {
0181 if (energy > threshold_) {
0182 int eta = cell.ieta();
0183 int phi = cell.iphi();
0184 double eta1 = (eta >= 0) ? (eta + 0.1) : (eta - 0.1);
0185 hdi2Oc_->Fill(eta1, (phi - 0.1));
0186 if (bh) {
0187 hdigEn_->Fill(energy);
0188 hdigOc_->Fill(eta1, (phi - 0.1));
0189 hdigLn_->Fill(depth);
0190 hdi3Oc_->Fill(eta1, depth);
0191 ++kount;
0192 edm::LogVerbatim("HGCalValidation")
0193 << "HGCalBHDigit[" << kount << "] ID " << cell << " E " << energy << ":" << (energy > threshold_);
0194 }
0195 }
0196 }
0197
0198 DEFINE_FWK_MODULE(HGCalBHValidation);
