HGCalValidation/test/HGCalRecHitStudy.cc

0001 // system include files
0002 #include <cmath>
0003 #include <iostream>
0004 #include <fstream>
0005 #include <map>
0006 #include <vector>
0007 #include <string>
0008
0009 // user include files
0010 #include "FWCore/Framework/interface/Frameworkfwd.h"
0011 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0012 #include "FWCore/Framework/interface/Event.h"
0013 #include "FWCore/Framework/interface/EventSetup.h"
0014 #include "FWCore/Framework/interface/ESHandle.h"
0015 #include "FWCore/Utilities/interface/InputTag.h"
0016 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0017 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0018 #include "FWCore/ServiceRegistry/interface/Service.h"
0019
0020 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0021
0022 #include "DataFormats/DetId/interface/DetId.h"
0023 #include "DataFormats/ForwardDetId/interface/HGCSiliconDetId.h"
0024 #include "DataFormats/ForwardDetId/interface/HGCScintillatorDetId.h"
0025 #include "DataFormats/HGCRecHit/interface/HGCRecHitCollections.h"
0026
0027 #include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
0028 #include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h"
0029 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0030
0031 #include "TH1D.h"
0032 #include "TH2D.h"
0033
0034 class HGCalRecHitStudy : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0035 public:
0036   explicit HGCalRecHitStudy(const edm::ParameterSet&);
0037   ~HGCalRecHitStudy() override = default;
0038
0039   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0040
0041 private:
0042   struct energysum {
0043     energysum() { e15 = e25 = e50 = e100 = e250 = e1000 = 0.0; }
0044     double e15, e25, e50, e100, e250, e1000;
0045   };
0046
0047   struct HitsInfo {
0048     HitsInfo() {
0049       x = y = z = time = energy = phi = eta = 0.0;
0050       layer = 0;
0051     }
0052     float x, y, z, time, energy, phi, eta;
0053     int layer;
0054   };
0055
0056   void beginJob() override {}
0057   void endJob() override {}
0058   void beginRun(edm::Run const&, edm::EventSetup const&) override;
0059   void endRun(edm::Run const&, edm::EventSetup const&) override {}
0060   void analyze(edm::Event const&, edm::EventSetup const&) override;
0061
0062   template <class T1, class T2>
0063   void recHitValidation(DetId& detId, int layer, const T1* geom, T2 it);
0064   void fillHitsInfo();
0065   void fillHitsInfo(HitsInfo& hits);
0066   void fillOccupancyMap(std::map<int, int>& occupancyMap, int layer);
0067
0068   // ----------member data ---------------------------
0069   const std::string nameDetector_;
0070   const edm::EDGetTokenT<HGCRecHitCollection> recHitSource_;
0071   const edm::ESGetToken<HGCalDDDConstants, IdealGeometryRecord> tok_hgcaldd_;
0072   const edm::ESGetToken<HGCalGeometry, IdealGeometryRecord> tok_hgcGeom_;
0073   const bool ifNose_, ifLayer_;
0074   const int verbosity_, nbinR_, nbinZ_, nbinEta_, nLayers_;
0075   const double rmin_, rmax_, zmin_, zmax_, etamin_, etamax_;
0076   unsigned int layers_;
0077   int firstLayer_;
0078   std::map<int, int> occupancyMapPlus_, occupancyMapMinus_;
0079
0080   std::vector<TH2D*> h_EtaPhiPlus_, h_EtaPhiMinus_, h_XY_;
0081   std::vector<TH1D*> h_energy_, h_HitOccupancyPlus_, h_HitOccupancyMinus_;
0082   TH1D *h_MeanHitOccupancyPlus_, *h_MeanHitOccupancyMinus_;
0083   TH2D *h_RZ_, *h_EtaPhi_;
0084 };
0085
0086 HGCalRecHitStudy::HGCalRecHitStudy(const edm::ParameterSet& iConfig)
0087     : nameDetector_(iConfig.getParameter<std::string>("detectorName")),
0088       recHitSource_(consumes<HGCRecHitCollection>(iConfig.getParameter<edm::InputTag>("source"))),
0089       tok_hgcaldd_(esConsumes<HGCalDDDConstants, IdealGeometryRecord, edm::Transition::BeginRun>(
0090           edm::ESInputTag{"", nameDetector_})),
0091       tok_hgcGeom_(esConsumes<HGCalGeometry, IdealGeometryRecord>(edm::ESInputTag{"", nameDetector_})),
0092       ifNose_(iConfig.getUntrackedParameter<bool>("ifNose", false)),
0093       ifLayer_(iConfig.getUntrackedParameter<bool>("ifLayer", false)),
0094       verbosity_(iConfig.getUntrackedParameter<int>("verbosity", 0)),
0095       nbinR_(iConfig.getUntrackedParameter<int>("nBinR", 300)),
0096       nbinZ_(iConfig.getUntrackedParameter<int>("nBinZ", 300)),
0097       nbinEta_(iConfig.getUntrackedParameter<int>("nBinEta", 100)),
0098       nLayers_(iConfig.getUntrackedParameter<int>("layers", 28)),
0099       rmin_(iConfig.getUntrackedParameter<double>("rMin", 0.0)),
0100       rmax_(iConfig.getUntrackedParameter<double>("rMax", 300.0)),
0101       zmin_(iConfig.getUntrackedParameter<double>("zMin", 300.0)),
0102       zmax_(iConfig.getUntrackedParameter<double>("zMax", 600.0)),
0103       etamin_(iConfig.getUntrackedParameter<double>("etaMin", 1.2)),
0104       etamax_(iConfig.getUntrackedParameter<double>("etaMax", 3.2)),
0105       layers_(0),
0106       firstLayer_(1) {
0107   usesResource(TFileService::kSharedResource);
0108
0109   edm::LogVerbatim("HGCalValidation") << "Initialize HGCalRecHitStudy for " << nameDetector_ << " with i/p tag "
0110                                       << iConfig.getParameter<edm::InputTag>("source") << " Flag " << ifNose_ << ":"
0111                                       << verbosity_;
0112 }
0113
0114 void HGCalRecHitStudy::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0115   edm::ParameterSetDescription desc;
0116   desc.add<std::string>("detectorName", "HGCalEESensitive");
0117   desc.add<edm::InputTag>("source", edm::InputTag("HGCalRecHit", "HGCEERecHits"));
0118   desc.addUntracked<bool>("ifNose", false);
0119   desc.addUntracked<bool>("ifLayer", false);
0120   desc.addUntracked<int>("verbosity", 0);
0121   desc.addUntracked<int>("nBinR", 300);
0122   desc.addUntracked<int>("nBinZ", 300);
0123   desc.addUntracked<int>("nBinEta", 200);
0124   desc.addUntracked<int>("layers", 28);
0125   desc.addUntracked<double>("rMin", 0.0);
0126   desc.addUntracked<double>("rMax", 300.0);
0127   desc.addUntracked<double>("zMin", 300.0);
0128   desc.addUntracked<double>("zMax", 600.0);
0129   desc.addUntracked<double>("etaMin", 1.0);
0130   desc.addUntracked<double>("etaMax", 3.0);
0131   descriptions.add("hgcalRecHitStudyEE", desc);
0132 }
0133
0134 void HGCalRecHitStudy::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0135   occupancyMapPlus_.clear();
0136   occupancyMapMinus_.clear();
0137
0138   bool ok(true);
0139   unsigned int ntot(0), nused(0);
0140   const edm::ESHandle<HGCalGeometry>& geom = iSetup.getHandle(tok_hgcGeom_);
0141   if (!geom.isValid())
0142     edm::LogWarning("HGCalValidation") << "Cannot get valid HGCalGeometry Object for " << nameDetector_;
0143   const HGCalGeometry* geom0 = geom.product();
0144
0145   const edm::Handle<HGCRecHitCollection>& theRecHitContainers = iEvent.getHandle(recHitSource_);
0146   if (theRecHitContainers.isValid()) {
0147     if (verbosity_ > 0)
0148       edm::LogVerbatim("HGCalValidation") << nameDetector_ << " with " << theRecHitContainers->size() << " element(s)";
0149     for (const auto& it : *(theRecHitContainers.product())) {
0150       ntot++;
0151       nused++;
0152       DetId detId = it.id();
0153       int layer = (ifNose_ ? HFNoseDetId(detId).layer()
0154                            : ((detId.det() == DetId::HGCalHSc) ? HGCScintillatorDetId(detId).layer()
0155                                                                : HGCSiliconDetId(detId).layer()));
0156       recHitValidation(detId, layer, geom0, &it);
0157     }
0158   } else {
0159     ok = false;
0160     edm::LogWarning("HGCalValidation") << "HGCRecHitCollection Handle does not exist !!!";
0161   }
0162   if (ok)
0163     fillHitsInfo();
0164   edm::LogVerbatim("HGCalValidation") << "Event " << iEvent.id().event() << " with " << ntot << " total and " << nused
0165                                       << " used recHits";
0166 }
0167
0168 template <class T1, class T2>
0169 void HGCalRecHitStudy::recHitValidation(DetId& detId, int layer, const T1* geom, T2 it) {
0170   GlobalPoint global = geom->getPosition(detId);
0171   double energy = it->energy();
0172
0173   float globalx = global.x();
0174   float globaly = global.y();
0175   float globalz = global.z();
0176   h_RZ_->Fill(std::abs(globalz), global.perp());
0177   if (ifLayer_) {
0178     if ((layer - firstLayer_) <= static_cast<int>(h_XY_.size()))
0179       h_XY_[layer - firstLayer_]->Fill(std::abs(globalx), std::abs(globaly));
0180   } else {
0181     h_EtaPhi_->Fill(std::abs(global.eta()), global.phi());
0182   }
0183   HitsInfo hinfo;
0184   hinfo.energy = energy;
0185   hinfo.x = globalx;
0186   hinfo.y = globaly;
0187   hinfo.z = globalz;
0188   hinfo.layer = layer - firstLayer_;
0189   hinfo.phi = global.phi();
0190   hinfo.eta = global.eta();
0191
0192   if (verbosity_ > 1)
0193     edm::LogVerbatim("HGCalValidation") << " --------------------------   gx = " << globalx << " gy = " << globaly
0194                                         << " gz = " << globalz << " phi = " << hinfo.phi << " eta = " << hinfo.eta;
0195
0196   fillHitsInfo(hinfo);
0197
0198   if (hinfo.eta > 0)
0199     fillOccupancyMap(occupancyMapPlus_, hinfo.layer);
0200   else
0201     fillOccupancyMap(occupancyMapMinus_, hinfo.layer);
0202 }
0203
0204 void HGCalRecHitStudy::fillOccupancyMap(std::map<int, int>& occupancyMap, int layer) {
0205   if (occupancyMap.find(layer) != occupancyMap.end())
0206     occupancyMap[layer]++;
0207   else
0208     occupancyMap[layer] = 1;
0209 }
0210
0211 void HGCalRecHitStudy::fillHitsInfo() {
0212   for (auto const& itr : occupancyMapPlus_) {
0213     int layer = itr.first;
0214     int occupancy = itr.second;
0215     h_HitOccupancyPlus_.at(layer)->Fill(occupancy);
0216   }
0217
0218   for (auto const& itr : occupancyMapMinus_) {
0219     int layer = itr.first;
0220     int occupancy = itr.second;
0221     h_HitOccupancyMinus_.at(layer)->Fill(occupancy);
0222   }
0223 }
0224
0225 void HGCalRecHitStudy::fillHitsInfo(HitsInfo& hits) {
0226   unsigned int ilayer = hits.layer;
0227   h_energy_.at(ilayer)->Fill(hits.energy);
0228   if (!ifLayer_) {
0229     h_EtaPhiPlus_.at(ilayer)->Fill(hits.eta, hits.phi);
0230     h_EtaPhiMinus_.at(ilayer)->Fill(hits.eta, hits.phi);
0231   }
0232 }
0233
0234 void HGCalRecHitStudy::beginRun(edm::Run const&, edm::EventSetup const& iSetup) {
0235   const HGCalDDDConstants& hgcons_ = iSetup.getData(tok_hgcaldd_);
0236   layers_ = hgcons_.layers(true);
0237   firstLayer_ = hgcons_.firstLayer();
0238
0239   edm::LogVerbatim("HGCalValidation") << "Finds " << layers_ << " layers for " << nameDetector_;
0240
0241   edm::Service<TFileService> fs;
0242   char histoname[100], title[100];
0243   for (unsigned int il = 0; il < layers_; il++) {
0244     int ilayer = firstLayer_ + static_cast<int>(il);
0245     sprintf(histoname, "HitOccupancy_Plus_layer_%d", ilayer);
0246     h_HitOccupancyPlus_.push_back(fs->make<TH1D>(histoname, "RecHitOccupancy_Plus", 100, 0, 10000));
0247     sprintf(histoname, "HitOccupancy_Minus_layer_%d", ilayer);
0248     h_HitOccupancyMinus_.push_back(fs->make<TH1D>(histoname, "RecHitOccupancy_Minus", 100, 0, 10000));
0249
0250     if (ifLayer_) {
0251       sprintf(histoname, "XY_L%d", ilayer);
0252       sprintf(title, "Y vs X at layer %d of %s", (ilayer + 1), nameDetector_.c_str());
0253       h_XY_.emplace_back(fs->make<TH2D>(histoname, title, nbinR_, 0, rmax_, nbinR_, 0, rmax_));
0254     } else {
0255       sprintf(histoname, "EtaPhi_Plus_Layer_%d", ilayer);
0256       h_EtaPhiPlus_.push_back(fs->make<TH2D>(histoname, "Occupancy", nbinEta_, etamin_, etamax_, 72, -M_PI, M_PI));
0257       sprintf(histoname, "EtaPhi_Minus_Layer_%d", ilayer);
0258       h_EtaPhiMinus_.push_back(fs->make<TH2D>(histoname, "Occupancy", nbinEta_, -etamax_, -etamin_, 72, -M_PI, M_PI));
0259     }
0260
0261     sprintf(histoname, "Energy_Layer_%d", ilayer);
0262     h_energy_.push_back(fs->make<TH1D>(histoname, "Energy", 1000, 0, 10.0));
0263   }  //loop over layers ends here
0264
0265   h_MeanHitOccupancyPlus_ =
0266       fs->make<TH1D>("SUMOfRecHitOccupancy_Plus", "SUMOfRecHitOccupancy_Plus", layers_, -0.5, layers_ - 0.5);
0267   h_MeanHitOccupancyMinus_ =
0268       fs->make<TH1D>("SUMOfRecHitOccupancy_Minus", "SUMOfRecHitOccupancy_Minus", layers_, -0.5, layers_ - 0.5);
0269
0270   sprintf(histoname, "RZ_%s", nameDetector_.c_str());
0271   sprintf(title, "R vs Z for %s", nameDetector_.c_str());
0272   h_RZ_ = fs->make<TH2D>(histoname, title, nbinZ_, zmin_, zmax_, nbinR_, rmin_, rmax_);
0273   if (!ifLayer_) {
0274     sprintf(histoname, "EtaPhi_%s", nameDetector_.c_str());
0275     sprintf(title, "#phi vs #eta for %s", nameDetector_.c_str());
0276     h_EtaPhi_ = fs->make<TH2D>(histoname, title, nbinEta_, etamin_, etamax_, 72, -M_PI, M_PI);
0277   }
0278 }
0279
0280 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0281
0282 #include "FWCore/Framework/interface/MakerMacros.h"
0283
0284 //define this as a plug-in
0285 DEFINE_FWK_MODULE(HGCalRecHitStudy);