HGCalBHValidation.cc

Go to the documentation of this file.

// system include files
#include <iostream>
#include <fstream>
#include <vector>
#include <map>
#include <string>
 
// user include files
#include "CommonTools/UtilAlgos/interface/TFileService.h"
#include "CondFormats/GeometryObjects/interface/HcalParameters.h"
#include "DataFormats/HcalDetId/interface/HcalDetId.h"
#include "DataFormats/HcalDigi/interface/HcalDigiCollections.h"
#include "DataFormats/HGCDigi/interface/HGCDigiCollections.h"
 
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/one/EDAnalyzer.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/ServiceRegistry/interface/Service.h"
 
#include "Geometry/Records/interface/HcalParametersRcd.h"
#include "SimDataFormats/CaloHit/interface/PCaloHit.h"
#include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
#include "DataFormats/HcalDetId/interface/HcalTestNumbering.h"
#include "DataFormats/ForwardDetId/interface/HGCScintillatorDetId.h"
 
// Root objects
#include "TROOT.h"
#include "TSystem.h"
#include "TFile.h"
#include "TH1.h"
#include "TH2.h"
 
class HGCalBHValidation : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
public:
  explicit HGCalBHValidation(const edm::ParameterSet& ps);
  ~HGCalBHValidation() override {}
 
  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
 
protected:
  void beginJob() override {}
  void beginRun(edm::Run const&, edm::EventSetup const&) override;
  void endRun(edm::Run const&, edm::EventSetup const&) override {}
  void analyze(edm::Event const&, edm::EventSetup const&) override;
  template <class T>
  void analyzeDigi(const T&, double const&, bool const&, int const&, unsigned int&);
 
private:
  edm::Service<TFileService> fs_;
  const std::string g4Label_, hcalHits_;
  const edm::InputTag hcalDigis_;
  const int iSample_, geomType_;
  const double threshold_;
  const bool ifHCAL_;
  edm::EDGetTokenT<edm::PCaloHitContainer> tok_hits_;
  edm::EDGetToken tok_hbhe_;
  int etaMax_;
 
  TH1D *hsimE1_, *hsimE2_, *hsimTm_;
  TH1D *hsimLn_, *hdigEn_, *hdigLn_;
  TH2D *hsimOc_, *hsi2Oc_, *hsi3Oc_;
  TH2D *hdigOc_, *hdi2Oc_, *hdi3Oc_;
};
 
HGCalBHValidation::HGCalBHValidation(const edm::ParameterSet& ps)
    : g4Label_(ps.getUntrackedParameter<std::string>("ModuleLabel", "g4SimHits")),
      hcalHits_((ps.getUntrackedParameter<std::string>("HitCollection", "HcalHits"))),
      hcalDigis_(ps.getUntrackedParameter<edm::InputTag>("DigiCollection")),
      iSample_(ps.getUntrackedParameter<int>("Sample", 5)),
      geomType_(ps.getUntrackedParameter<int>("GeometryType", 0)),
      threshold_(ps.getUntrackedParameter<double>("Threshold", 12.0)),
      ifHCAL_(ps.getUntrackedParameter<bool>("ifHCAL", false)),
      etaMax_(100) {
  usesResource(TFileService::kSharedResource);
 
  tok_hits_ = consumes<edm::PCaloHitContainer>(edm::InputTag(g4Label_, hcalHits_));
  if (ifHCAL_)
    tok_hbhe_ = consumes<QIE11DigiCollection>(hcalDigis_);
  else
    tok_hbhe_ = consumes<HGCalDigiCollection>(hcalDigis_);
  edm::LogVerbatim("HGCalValidation") << "HGCalBHValidation::Input for SimHit:" << edm::InputTag(g4Label_, hcalHits_)
                                      << "  Digits:" << hcalDigis_ << "  Sample: " << iSample_ << "  Threshold "
                                      << threshold_;
}
 
void HGCalBHValidation::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
  edm::ParameterSetDescription desc;
  desc.addUntracked<std::string>("ModuleLabel", "g4SimHits");
  desc.addUntracked<std::string>("HitCollection", "HGCHitsHEback");
  desc.addUntracked<edm::InputTag>("DigiCollection", edm::InputTag("simHGCalUnsuppressedDigis", "HEback"));
  desc.addUntracked<int>("Sample", 5);
  desc.addUntracked<int>("GeometryType", 1);
  desc.addUntracked<double>("Threshold", 15.0);
  desc.addUntracked<bool>("ifHCAL", false);
  descriptions.add("hgcalBHAnalysis", desc);
}
 
void HGCalBHValidation::beginRun(edm::Run const&, edm::EventSetup const& es) {
  if (geomType_ == 0) {
    std::string label;
    edm::ESHandle<HcalParameters> parHandle;
    es.get<HcalParametersRcd>().get(label, parHandle);
    const HcalParameters* hpar = &(*parHandle);
    const std::vector<int> etaM = hpar->etaMax;
    etaMax_ = etaM[1];
  }
  edm::LogVerbatim("HGCalValidation") << "HGCalBHValidation::Maximum Number of"
                                      << " eta sectors:" << etaMax_ << "\nHitsValidationHcal::Booking the Histograms";
 
  //Histograms for Sim Hits
  hsimE1_ = fs_->make<TH1D>("SimHitEn1", "Sim Hit Energy", 1000, 0.0, 1.0);
  hsimE2_ = fs_->make<TH1D>("SimHitEn2", "Sim Hit Energy", 1000, 0.0, 1.0);
  hsimTm_ = fs_->make<TH1D>("SimHitTime", "Sim Hit Time", 1000, 0.0, 500.0);
  hsimLn_ = fs_->make<TH1D>("SimHitLong", "Sim Hit Long. Profile", 50, 0.0, 25.0);
  hsimOc_ = fs_->make<TH2D>("SimHitOccup", "Sim Hit Occupnacy", 2 * etaMax_ + 1, -etaMax_, etaMax_ + 1, 360, 0, 360);
  hsi2Oc_ = fs_->make<TH2D>("SimHitOccu2", "Sim Hit Occupnacy", 2 * etaMax_ + 1, -etaMax_, etaMax_ + 1, 360, 0, 360);
  hsi3Oc_ = fs_->make<TH2D>("SimHitOccu3", "Sim Hit Occupnacy", 2 * etaMax_ + 1, -etaMax_, etaMax_ + 1, 50, 0, 25);
  //Histograms for Digis
  hdigEn_ = fs_->make<TH1D>("DigiEnergy", "Digi ADC Sample", 1000, 0.0, 1000.0);
  hdigLn_ = fs_->make<TH1D>("DigiLong", "Digi Long. Profile", 50, 0.0, 25.0);
  hdigOc_ = fs_->make<TH2D>("DigiOccup", "Digi Occupnacy", 2 * etaMax_ + 1, -etaMax_, etaMax_ + 1, 360, 0, 360);
  hdi2Oc_ = fs_->make<TH2D>("DigiOccu2", "Digi Occupnacy", 2 * etaMax_ + 1, -etaMax_, etaMax_ + 1, 360, 0, 360);
  hdi3Oc_ = fs_->make<TH2D>("DigiOccu3", "Digi Occupnacy", 2 * etaMax_ + 1, -etaMax_, etaMax_ + 1, 50, 0, 25);
}
 
void HGCalBHValidation::analyze(const edm::Event& e, const edm::EventSetup&) {
  edm::LogVerbatim("HGCalValidation") << "HGCalBHValidation:Run = " << e.id().run() << " Event = " << e.id().event();
 
  //SimHits
  edm::Handle<edm::PCaloHitContainer> hitsHcal;
  e.getByToken(tok_hits_, hitsHcal);
  edm::LogVerbatim("HGCalValidation") << "HGCalBHValidation.: PCaloHitContainer"
                                      << " obtained with flag " << hitsHcal.isValid();
  if (hitsHcal.isValid()) {
    edm::LogVerbatim("HGCalValidation") << "HGCalBHValidation: PCaloHit buffer " << hitsHcal->size();
    unsigned i(0);
    std::map<unsigned int, double> map_try;
    for (edm::PCaloHitContainer::const_iterator it = hitsHcal->begin(); it != hitsHcal->end(); ++it) {
      double energy = it->energy();
      double time = it->time();
      unsigned int id = it->id();
      int subdet(0), z(0), depth(0), eta(0), phi(0), lay(0);
      bool hbhe(false), bh(false);
      if (geomType_ == 0) {
        HcalTestNumbering::unpackHcalIndex(id, subdet, z, depth, eta, phi, lay);
        if (z == 0)
          eta = -eta;
        hbhe = ((subdet == static_cast<int>(HcalEndcap)) || (subdet == static_cast<int>(HcalBarrel)));
        bh = (subdet == static_cast<int>(HcalEndcap));
      } else {
        hbhe = bh = (DetId(id).det() == DetId::HGCalHSc);
        if (bh) {
          eta = HGCScintillatorDetId(id).ieta();
          phi = HGCScintillatorDetId(id).iphi();
          lay = HGCScintillatorDetId(id).layer();
        }
      }
      double eta1 = (eta >= 0) ? (eta + 0.1) : (eta - 0.1);
      if (hbhe)
        hsi2Oc_->Fill(eta1, (phi - 0.1), energy);
      if (bh) {
        hsimE1_->Fill(energy);
        hsimTm_->Fill(time, energy);
        hsimOc_->Fill(eta1, (phi - 0.1), energy);
        hsimLn_->Fill(lay, energy);
        hsi3Oc_->Fill(eta1, lay, energy);
        double ensum(0);
        if (map_try.count(id) != 0)
          ensum = map_try[id];
        ensum += energy;
        map_try[id] = ensum;
        ++i;
        edm::LogVerbatim("HGCalValidation")
            << "HGCalBHHit[" << i << "] ID " << std::hex << " " << id << std::dec << " SubDet " << subdet << " depth "
            << depth << " Eta " << eta << " Phi " << phi << " layer " << lay << " E " << energy << " time " << time;
      }
    }
    for (std::map<unsigned int, double>::iterator itr = map_try.begin(); itr != map_try.end(); ++itr) {
      hsimE2_->Fill((*itr).second);
    }
  }
 
  //Digits
  unsigned int kount(0);
  if ((geomType_ == 0) && ifHCAL_) {
    edm::Handle<QIE11DigiCollection> hbhecoll;
    e.getByToken(tok_hbhe_, hbhecoll);
    edm::LogVerbatim("HGCalValidation") << "HGCalBHValidation.: "
                                        << "HBHEQIE11DigiCollection obtained "
                                        << "with flag " << hbhecoll.isValid();
    if (hbhecoll.isValid()) {
      edm::LogVerbatim("HGCalValidation") << "HGCalBHValidation: HBHEDigit "
                                          << "buffer " << hbhecoll->size();
      for (QIE11DigiCollection::const_iterator it = hbhecoll->begin(); it != hbhecoll->end(); ++it) {
        QIE11DataFrame df(*it);
        HcalDetId cell(df.id());
        bool hbhe =
            ((cell.subdetId() == static_cast<int>(HcalEndcap)) || (cell.subdetId() == static_cast<int>(HcalBarrel)));
        if (hbhe) {
          bool bh = (cell.subdetId() == static_cast<int>(HcalEndcap));
          int depth = cell.depth();
          double energy = df[iSample_].adc();
          analyzeDigi(cell, energy, bh, depth, kount);
        }
      }
    }
  } else {
    edm::Handle<HGCalDigiCollection> hbhecoll;
    e.getByToken(tok_hbhe_, hbhecoll);
    edm::LogVerbatim("HGCalValidation") << "HGCalBHValidation.: "
                                        << "HGCalDigiCollection obtained with"
                                        << " flag " << hbhecoll.isValid();
    if (hbhecoll.isValid()) {
      edm::LogVerbatim("HGCalValidation") << "HGCalBHValidation: HGCalDigi "
                                          << "buffer " << hbhecoll->size();
      for (HGCalDigiCollection::const_iterator it = hbhecoll->begin(); it != hbhecoll->end(); ++it) {
        HGCalDataFrame df(*it);
        double energy = df[iSample_].data();
        if (geomType_ == 0) {
          HcalDetId cell(df.id());
          bool hbhe =
              ((cell.subdetId() == static_cast<int>(HcalEndcap)) || (cell.subdetId() == static_cast<int>(HcalBarrel)));
          if (hbhe) {
            bool bh = (cell.subdetId() == static_cast<int>(HcalEndcap));
            int depth = cell.depth();
            analyzeDigi(cell, energy, bh, depth, kount);
          }
        } else {
          bool bh = (DetId(df.id()).det() == DetId::HGCalHSc);
          if (bh) {
            HGCScintillatorDetId cell(df.id());
            int depth = cell.layer();
            analyzeDigi(cell, energy, bh, depth, kount);
          }
        }
      }
    }
  }
}
 
template <class T>
void HGCalBHValidation::analyzeDigi(
    const T& cell, double const& energy, bool const& bh, int const& depth, unsigned int& kount) {
  if (energy > threshold_) {
    int eta = cell.ieta();
    int phi = cell.iphi();
    double eta1 = (eta >= 0) ? (eta + 0.1) : (eta - 0.1);
    hdi2Oc_->Fill(eta1, (phi - 0.1));
    if (bh) {
      hdigEn_->Fill(energy);
      hdigOc_->Fill(eta1, (phi - 0.1));
      hdigLn_->Fill(depth);
      hdi3Oc_->Fill(eta1, depth);
      ++kount;
      edm::LogVerbatim("HGCalValidation")
          << "HGCalBHDigit[" << kount << "] ID " << cell << " E " << energy << ":" << (energy > threshold_);
    }
  }
}
 
DEFINE_FWK_MODULE(HGCalBHValidation);