Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-05-05 03:16:18

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 {}
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   //Histograms for Sim Hits
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   //Histograms for Digis
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   //SimHits
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   //Digits
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);