Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-09-21 22:40:25

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/ESHandle.h"
0017 #include "FWCore/Framework/interface/MakerMacros.h"
0018 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0019 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0020 #include "FWCore/ServiceRegistry/interface/Service.h"
0021 
0022 #include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
0023 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0024 
0025 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0026 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0027 #include "DataFormats/ForwardDetId/interface/HGCSiliconDetId.h"
0028 
0029 // Root objects
0030 #include "TROOT.h"
0031 #include "TSystem.h"
0032 #include "TFile.h"
0033 #include "TH1.h"
0034 #include "TH2.h"
0035 
0036 class HGCalSiliconValidation : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0037 public:
0038   explicit HGCalSiliconValidation(const edm::ParameterSet& ps);
0039   ~HGCalSiliconValidation() override = default;
0040 
0041   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0042 
0043 protected:
0044   void beginJob() override {}
0045   void beginRun(edm::Run const&, edm::EventSetup const&) override;
0046   void endRun(edm::Run const&, edm::EventSetup const&) override {}
0047   void analyze(edm::Event const&, edm::EventSetup const&) override;
0048 
0049 private:
0050   const std::string g4Label_, nameDetector_, hgcalHits_;
0051   const edm::InputTag hgcalDigis_;
0052   const int iSample_;
0053   const edm::ESGetToken<HGCalGeometry, IdealGeometryRecord> tok_hgcalgeom_;
0054   const edm::EDGetTokenT<edm::PCaloHitContainer> tok_hits_;
0055   const edm::EDGetTokenT<HGCalDigiCollection> tok_digi_;
0056 
0057   TH1D *hsimE1_, *hsimE2_, *hsimTm_;
0058   TH1D *hsimLn_, *hdigEn_, *hdigLn_;
0059   TH2D *hsimOc_, *hsi2Oc_, *hdigOc_, *hdi2Oc_;
0060 };
0061 
0062 HGCalSiliconValidation::HGCalSiliconValidation(const edm::ParameterSet& ps)
0063     : g4Label_(ps.getUntrackedParameter<std::string>("ModuleLabel", "g4SimHits")),
0064       nameDetector_(ps.getUntrackedParameter<std::string>("detectorName", "HGCalEESensitive")),
0065       hgcalHits_((ps.getUntrackedParameter<std::string>("HitCollection", "HGCHitsEE"))),
0066       hgcalDigis_(ps.getUntrackedParameter<edm::InputTag>("DigiCollection")),
0067       iSample_(ps.getUntrackedParameter<int>("Sample", 5)),
0068       tok_hgcalgeom_(esConsumes<HGCalGeometry, IdealGeometryRecord>(edm::ESInputTag{"", nameDetector_})),
0069       tok_hits_(consumes<edm::PCaloHitContainer>(edm::InputTag(g4Label_, hgcalHits_))),
0070       tok_digi_(consumes<HGCalDigiCollection>(hgcalDigis_)) {
0071   usesResource(TFileService::kSharedResource);
0072 
0073   edm::LogVerbatim("HGCalValidation") << "HGCalSiliconValidation::Input for SimHit:"
0074                                       << edm::InputTag(g4Label_, hgcalHits_) << "  Digits:" << hgcalDigis_
0075                                       << "  Sample: " << iSample_;
0076 }
0077 
0078 void HGCalSiliconValidation::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0079   edm::ParameterSetDescription desc;
0080   desc.addUntracked<std::string>("ModuleLabel", "g4SimHits");
0081   desc.addUntracked<std::string>("detectorName", "HGCalEESensitive");
0082   desc.addUntracked<std::string>("HitCollection", "HGCHitsEE");
0083   desc.addUntracked<edm::InputTag>("DigiCollection", edm::InputTag("simHGCalUnsuppressedDigis", "EE"));
0084   desc.addUntracked<int>("Sample", 5);
0085   descriptions.add("hgcalSiliconAnalysisEE", desc);
0086 }
0087 
0088 void HGCalSiliconValidation::beginRun(edm::Run const&, edm::EventSetup const& es) {
0089   edm::LogVerbatim("HGCalValidation") << "HGCalSiliconValidation::Booking the Histograms";
0090 
0091   //Histograms for Sim Hits
0092   edm::Service<TFileService> fs;
0093   hsimE1_ = fs->make<TH1D>("SimHitEn1", "Sim Hit Energy", 1000, 0.0, 1.0);
0094   hsimE2_ = fs->make<TH1D>("SimHitEn2", "Sim Hit Energy", 1000, 0.0, 1.0);
0095   hsimTm_ = fs->make<TH1D>("SimHitTime", "Sim Hit Time", 1000, 0.0, 500.0);
0096   hsimLn_ = fs->make<TH1D>("SimHitLong", "Sim Hit Long. Profile", 60, 0.0, 30.0);
0097   hsimOc_ = fs->make<TH2D>("SimHitOccup", "Sim Hit Occupnacy", 300, 0.0, 300.0, 60, 0.0, 30.0);
0098   hsi2Oc_ = fs->make<TH2D>("SimHitOccu2", "Sim Hit Occupnacy", 300, 300.0, 600.0, 300, 0.0, 300.0);
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", 60, 0.0, 30.0);
0102   hdigOc_ = fs->make<TH2D>("DigiOccup", "Digi Occupnacy", 300, 0.0, 300.0, 60, 0.0, 30.0);
0103   hdi2Oc_ = fs->make<TH2D>("DigiOccu2", "Digi Occupnacy", 300, 300.0, 600.0, 300, 0.0, 300.0);
0104 }
0105 
0106 void HGCalSiliconValidation::analyze(const edm::Event& e, const edm::EventSetup& iSetup) {
0107   edm::LogVerbatim("HGCalValidation") << "HGCalSiliconValidation:Run = " << e.id().run()
0108                                       << " Event = " << e.id().event();
0109 
0110   const edm::ESHandle<HGCalGeometry>& geom = iSetup.getHandle(tok_hgcalgeom_);
0111   if (!geom.isValid()) {
0112     edm::LogWarning("HGCalValidation") << "Cannot get valid HGCalGeometry Object for " << nameDetector_;
0113   } else {
0114     const HGCalGeometry* geom0 = geom.product();
0115 
0116     //SimHits
0117     const edm::Handle<edm::PCaloHitContainer>& hitsCalo = e.getHandle(tok_hits_);
0118     edm::LogVerbatim("HGCalValidation") << "HGCalSiliconValidation.: PCaloHitContainer obtained with flag "
0119                                         << hitsCalo.isValid();
0120     if (hitsCalo.isValid()) {
0121       edm::LogVerbatim("HGCalValidation") << "HGCalSiliconValidation: PCaloHit buffer " << hitsCalo->size();
0122       unsigned i(0);
0123       std::map<unsigned int, double> map_try;
0124       for (edm::PCaloHitContainer::const_iterator it = hitsCalo->begin(); it != hitsCalo->end(); ++it) {
0125         double energy = it->energy();
0126         double time = it->time();
0127         unsigned int id = it->id();
0128         GlobalPoint pos = geom0->getPosition(DetId(id));
0129         double r = pos.perp();
0130         double z = std::abs(pos.z());
0131         int lay = HGCSiliconDetId(id).layer();
0132         hsimE1_->Fill(energy);
0133         hsimTm_->Fill(time, energy);
0134         hsimOc_->Fill(r, lay, energy);
0135         hsi2Oc_->Fill(z, r, energy);
0136         hsimLn_->Fill(lay, energy);
0137         double ensum = (map_try.count(id) != 0) ? map_try[id] : 0;
0138         ensum += energy;
0139         map_try[id] = ensum;
0140         ++i;
0141         edm::LogVerbatim("HGCalValidation") << "HGCalBHHit[" << i << "] ID " << std::hex << " " << id << std::dec << " "
0142                                             << HGCSiliconDetId(id) << " E " << energy << " time " << time;
0143       }
0144       for (std::map<unsigned int, double>::iterator itr = map_try.begin(); itr != map_try.end(); ++itr) {
0145         hsimE2_->Fill((*itr).second);
0146       }
0147     }
0148 
0149     //Digits
0150     unsigned int kount(0);
0151     const edm::Handle<HGCalDigiCollection>& digicoll = e.getHandle(tok_digi_);
0152     edm::LogVerbatim("HGCalValidation") << "HGCalSiliconValidation.: HGCalDigiCollection obtained with flag "
0153                                         << digicoll.isValid();
0154     if (digicoll.isValid()) {
0155       edm::LogVerbatim("HGCalValidation") << "HGCalSiliconValidation: HGCalDigi buffer " << digicoll->size();
0156       for (HGCalDigiCollection::const_iterator it = digicoll->begin(); it != digicoll->end(); ++it) {
0157         HGCalDataFrame df(*it);
0158         double energy = df[iSample_].data();
0159         HGCSiliconDetId cell(df.id());
0160         GlobalPoint pos = geom0->getPosition(cell);
0161         double r = pos.perp();
0162         double z = std::abs(pos.z());
0163         int depth = cell.layer();
0164         hdigEn_->Fill(energy);
0165         hdigLn_->Fill(depth);
0166         hdigOc_->Fill(r, depth);
0167         hdi2Oc_->Fill(z, r);
0168         ++kount;
0169         edm::LogVerbatim("HGCalValidation") << "HGCalBHDigit[" << kount << "] ID " << cell << " E " << energy;
0170       }
0171     }
0172   }
0173 }
0174 
0175 DEFINE_FWK_MODULE(HGCalSiliconValidation);