Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-12-03 01:10:33

0001 // system include files
0002 #include <iostream>
0003 #include <fstream>
0004 #include <vector>
0005 #include <map>
0006 #include <string>
0007 
0008 // user include files
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 // Root objects
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   //Histograms for Sim Hits
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   //Histograms for Digis
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   //SimHits
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   //Digits
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);