HGCalTestBeam/plugins/HGCalTimingAnalyzer.cc

0001 // system include files
0002 #include <fstream>
0003 #include <sstream>
0004 #include <map>
0005 #include <string>
0006 #include <vector>
0007
0008 // user include files
0009 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0010 #include "DataFormats/ForwardDetId/interface/HGCalDetId.h"
0011 #include "DataFormats/HGCDigi/interface/HGCDigiCollections.h"
0012 #include "DataFormats/HGCRecHit/interface/HGCRecHitCollections.h"
0013
0014 #include "FWCore/Framework/interface/Event.h"
0015 #include "FWCore/Framework/interface/EventSetup.h"
0016 #include "FWCore/Framework/interface/Frameworkfwd.h"
0017 #include "FWCore/Framework/interface/MakerMacros.h"
0018 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0019 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0020 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0021 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0022 #include "FWCore/ServiceRegistry/interface/Service.h"
0023 #include "FWCore/Utilities/interface/InputTag.h"
0024
0025 #include "Geometry/HGCalTBCommonData/interface/HGCalTBDDDConstants.h"
0026 #include "Geometry/HGCalGeometry/interface/HGCalTBGeometry.h"
0027 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0028 #include "Geometry/Records/interface/HcalRecNumberingRecord.h"
0029 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0030 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0031 #include "SimDataFormats/CaloTest/interface/HGCalTestNumbering.h"
0032 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0033 #include "SimDataFormats/HcalTestBeam/interface/HcalTestBeamNumbering.h"
0034 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0035 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
0036
0037 // Root objects
0038 #include "TROOT.h"
0039 #include "TSystem.h"
0040 #include "TTree.h"
0041
0042 //#define EDM_ML_DEBUG
0043
0044 class HGCalTimingAnalyzer : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0045 public:
0046   explicit HGCalTimingAnalyzer(edm::ParameterSet const&);
0047   ~HGCalTimingAnalyzer() override = default;
0048
0049   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0050
0051 private:
0052   void beginJob() override;
0053   void beginRun(edm::Run const&, edm::EventSetup const&) override;
0054   void endRun(edm::Run const&, edm::EventSetup const&) override {}
0055   void analyze(edm::Event const&, edm::EventSetup const&) override;
0056   void analyzeSimHits(int type, std::vector<PCaloHit> const& hits);
0057   void analyzeSimTracks(edm::Handle<edm::SimTrackContainer> const& SimTk,
0058                         edm::Handle<edm::SimVertexContainer> const& SimVtx);
0059
0060   edm::Service<TFileService> fs_;
0061   const std::vector<int> idBeamDef_ = {1001};
0062   const std::string detectorEE_, detectorBeam_;
0063   const bool groupHits_;
0064   const double timeUnit_;
0065   const bool doTree_;
0066   const std::vector<int> idBeams_;
0067   const edm::ESGetToken<HGCalTBDDDConstants, IdealGeometryRecord> tokDDD_;
0068   const HGCalTBDDDConstants* hgcons_;
0069   const edm::InputTag labelGen_;
0070   const std::string labelHitEE_, labelHitBeam_;
0071   const edm::EDGetTokenT<edm::HepMCProduct> tok_hepMC_;
0072   const edm::EDGetTokenT<edm::SimTrackContainer> tok_simTk_;
0073   const edm::EDGetTokenT<edm::SimVertexContainer> tok_simVtx_;
0074   const edm::EDGetTokenT<edm::PCaloHitContainer> tok_hitsEE_, tok_hitsBeam_;
0075   TTree* tree_;
0076   std::vector<uint32_t> simHitCellIdEE_, simHitCellIdBeam_;
0077   std::vector<float> simHitCellEnEE_, simHitCellEnBeam_;
0078   std::vector<float> simHitCellTmEE_, simHitCellTmBeam_;
0079   double xBeam_, yBeam_, zBeam_, pBeam_;
0080 };
0081
0082 HGCalTimingAnalyzer::HGCalTimingAnalyzer(const edm::ParameterSet& iConfig)
0083     : detectorEE_(iConfig.getParameter<std::string>("DetectorEE")),
0084       detectorBeam_(iConfig.getParameter<std::string>("DetectorBeam")),
0085       groupHits_(iConfig.getParameter<bool>("GroupHits")),
0086       timeUnit_((!groupHits_) ? 0.000001 : (iConfig.getParameter<double>("TimeUnit"))),
0087       doTree_(iConfig.getUntrackedParameter<bool>("DoTree", false)),
0088       idBeams_((iConfig.getParameter<std::vector<int>>("IDBeams")).empty()
0089                    ? idBeamDef_
0090                    : (iConfig.getParameter<std::vector<int>>("IDBeams"))),
0091       tokDDD_(esConsumes<HGCalTBDDDConstants, IdealGeometryRecord, edm::Transition::BeginRun>(
0092           edm::ESInputTag("", detectorEE_))),
0093       labelGen_(iConfig.getParameter<edm::InputTag>("GeneratorSrc")),
0094       labelHitEE_(iConfig.getParameter<std::string>("CaloHitSrcEE")),
0095       labelHitBeam_(iConfig.getParameter<std::string>("CaloHitSrcBeam")),
0096       tok_hepMC_(consumes<edm::HepMCProduct>(labelGen_)),
0097       tok_simTk_(consumes<edm::SimTrackContainer>(edm::InputTag("g4SimHits"))),
0098       tok_simVtx_(consumes<edm::SimVertexContainer>(edm::InputTag("g4SimHits"))),
0099       tok_hitsEE_(consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits", labelHitEE_))),
0100       tok_hitsBeam_(consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits", labelHitBeam_))) {
0101   usesResource("TFileService");
0102
0103   // now do whatever initialization is needed
0104   // Group hits (if groupHits_ = true) if hits come within timeUnit_
0105   // Only look into the beam counters with ID's as in idBeams_
0106 #ifdef EDM_ML_DEBUG
0107   std::ostringstream st1;
0108   st1 << "HGCalTimingAnalyzer:: Group Hits " << groupHits_ << " in " << timeUnit_ << " IdBeam " << idBeams_.size()
0109       << ":";
0110   for (const auto& id : idBeams_)
0111     st1 << " " << id;
0112   edm::LogVerbatim("HGCSim") << st1.str();
0113
0114   edm::LogVerbatim("HGCSim") << "HGCalTimingAnalyzer:: GeneratorSource = " << labelGen_;
0115   edm::LogVerbatim("HGCSim") << "HGCalTimingAnalyzer:: Detector " << detectorEE_ << " with tags " << labelHitEE_;
0116   edm::LogVerbatim("HGCSim") << "HGCalTimingAnalyzer:: Detector " << detectorBeam_ << " with tags " << labelHitBeam_;
0117 #endif
0118 }
0119
0120 void HGCalTimingAnalyzer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0121   edm::ParameterSetDescription desc;
0122   desc.add<std::string>("DetectorEE", "HGCalEESensitive");
0123   desc.add<std::string>("DetectorBeam", "HcalTB06BeamDetector");
0124   desc.add<bool>("GroupHits", false);
0125   desc.add<double>("TimeUnit", 0.001);
0126   std::vector<int> ids = {1001, 1002, 1003, 1004, 1005};
0127   desc.add<std::vector<int>>("IDBeams", ids);
0128   desc.addUntracked<bool>("DoTree", true);
0129   desc.add<edm::InputTag>("GeneratorSrc", edm::InputTag("generatorSmeared"));
0130   desc.add<std::string>("CaloHitSrcEE", "HGCHitsEE");
0131   desc.add<std::string>("CaloHitSrcBeam", "HcalTB06BeamHits");
0132   descriptions.add("HGCalTimingAnalyzer", desc);
0133 }
0134
0135 void HGCalTimingAnalyzer::beginJob() {
0136   std::string det(detectorEE_);
0137   if (doTree_) {
0138     tree_ = fs_->make<TTree>("HGCTB", "SimHitEnergy");
0139     tree_->Branch("xBeam", &xBeam_, "xBeam/D");
0140     tree_->Branch("yBeam", &yBeam_, "yBeam/D");
0141     tree_->Branch("zBeam", &zBeam_, "zBeam/D");
0142     tree_->Branch("pBeam", &pBeam_, "pBeam/D");
0143     tree_->Branch("simHitCellIdEE_", &simHitCellIdEE_);
0144     tree_->Branch("simHitCellEnEE_", &simHitCellEnEE_);
0145     tree_->Branch("simHitCellTmEE_", &simHitCellTmEE_);
0146     tree_->Branch("simHitCellIdBeam_", &simHitCellIdBeam_);
0147     tree_->Branch("simHitCellEnBeam_", &simHitCellEnBeam_);
0148     tree_->Branch("simHitCellTmBeam_", &simHitCellTmBeam_);
0149   }
0150 }
0151
0152 void HGCalTimingAnalyzer::beginRun(const edm::Run&, const edm::EventSetup& iSetup) {
0153   hgcons_ = &iSetup.getData(tokDDD_);
0154
0155 #ifdef EDM_ML_DEBUG
0156   edm::LogVerbatim("HGCSim") << "HGCalTimingAnalyzer::" << detectorEE_ << " defined with " << hgcons_->layers(false)
0157                              << " layers";
0158 #endif
0159 }
0160
0161 void HGCalTimingAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0162 #ifdef EDM_ML_DEBUG
0163   // Generator input
0164   const edm::Handle<edm::HepMCProduct>& evtMC = iEvent.getHandle(tok_hepMC_);
0165   if (!evtMC.isValid()) {
0166     edm::LogWarning("HGCal") << "no HepMCProduct found";
0167   } else {
0168     const HepMC::GenEvent* myGenEvent = evtMC->GetEvent();
0169     unsigned int k(0);
0170     for (HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end();
0171          ++p, ++k) {
0172       edm::LogVerbatim("HGCSim") << "Particle[" << k << "] with p " << (*p)->momentum().rho() << " theta "
0173                                  << (*p)->momentum().theta() << " phi " << (*p)->momentum().phi();
0174     }
0175   }
0176 #endif
0177
0178   // Now the Simhits
0179   const edm::Handle<edm::SimTrackContainer>& SimTk = iEvent.getHandle(tok_simTk_);
0180   const edm::Handle<edm::SimVertexContainer>& SimVtx = iEvent.getHandle(tok_simVtx_);
0181   analyzeSimTracks(SimTk, SimVtx);
0182
0183   simHitCellIdEE_.clear();
0184   simHitCellIdBeam_.clear();
0185   simHitCellEnEE_.clear();
0186   simHitCellEnBeam_.clear();
0187   simHitCellTmEE_.clear();
0188   simHitCellTmBeam_.clear();
0189
0190   std::vector<PCaloHit> caloHits;
0191   const edm::Handle<edm::PCaloHitContainer>& theCaloHitContainers = iEvent.getHandle(tok_hitsEE_);
0192   if (theCaloHitContainers.isValid()) {
0193 #ifdef EDM_ML_DEBUG
0194     edm::LogVerbatim("HGCSim") << "PcalohitContainer for " << detectorEE_ << " has " << theCaloHitContainers->size()
0195                                << " hits";
0196 #endif
0197     caloHits.clear();
0198     caloHits.insert(caloHits.end(), theCaloHitContainers->begin(), theCaloHitContainers->end());
0199     analyzeSimHits(0, caloHits);
0200   } else {
0201 #ifdef EDM_ML_DEBUG
0202     edm::LogVerbatim("HGCSim") << "PCaloHitContainer does not exist for " << detectorEE_ << " !!!";
0203 #endif
0204   }
0205
0206   const edm::Handle<edm::PCaloHitContainer>& caloHitContainerBeam = iEvent.getHandle(tok_hitsBeam_);
0207   if (caloHitContainerBeam.isValid()) {
0208 #ifdef EDM_ML_DEBUG
0209     edm::LogVerbatim("HGCSim") << "PcalohitContainer for " << detectorBeam_ << " has " << caloHitContainerBeam->size()
0210                                << " hits";
0211 #endif
0212     caloHits.clear();
0213     caloHits.insert(caloHits.end(), caloHitContainerBeam->begin(), caloHitContainerBeam->end());
0214     analyzeSimHits(1, caloHits);
0215   } else {
0216 #ifdef EDM_ML_DEBUG
0217     edm::LogVerbatim("HGCSim") << "PCaloHitContainer does not exist for " << detectorBeam_ << " !!!";
0218 #endif
0219   }
0220   if (doTree_)
0221     tree_->Fill();
0222 }
0223
0224 void HGCalTimingAnalyzer::analyzeSimHits(int type, std::vector<PCaloHit> const& hits) {
0225 #ifdef EDM_ML_DEBUG
0226   unsigned int i(0);
0227 #endif
0228   std::map<std::pair<uint32_t, uint64_t>, std::pair<double, double>> map_hits;
0229   for (const auto& hit : hits) {
0230     double energy = hit.energy();
0231     double time = hit.time();
0232     uint32_t id = hit.id();
0233     if (type == 0) {
0234       int subdet, zside, layer, sector, subsector, cell;
0235       HGCalTestNumbering::unpackHexagonIndex(id, subdet, zside, layer, sector, subsector, cell);
0236       std::pair<int, int> recoLayerCell = hgcons_->simToReco(cell, layer, sector, true);
0237       id = HGCalDetId((ForwardSubdetector)(subdet), zside, recoLayerCell.second, subsector, sector, recoLayerCell.first)
0238                .rawId();
0239 #ifdef EDM_ML_DEBUG
0240       edm::LogVerbatim("HGCSim") << "SimHit:Hit[" << i << "] Id " << subdet << ":" << zside << ":" << layer << ":"
0241                                  << sector << ":" << subsector << ":" << recoLayerCell.first << ":"
0242                                  << recoLayerCell.second << " Energy " << energy << " Time " << time;
0243 #endif
0244     } else {
0245 #ifdef EDM_ML_DEBUG
0246       int subdet, layer, x, y;
0247       HcalTestBeamNumbering::unpackIndex(id, subdet, layer, x, y);
0248       edm::LogVerbatim("HGCSim") << "SimHit:Hit[" << i << "] Beam Subdet " << subdet << " Layer " << layer << " x|y "
0249                                  << x << ":" << y << " Energy " << energy << " Time " << time;
0250 #endif
0251     }
0252     uint64_t tid = (uint64_t)((time + 50.0) / timeUnit_);
0253     std::pair<uint32_t, uint64_t> key(id, tid);
0254     auto itr = map_hits.find(key);
0255     if (itr == map_hits.end()) {
0256       map_hits[key] = std::pair<double, double>(time, 0.0);
0257       itr = map_hits.find(key);
0258     }
0259     energy += (itr->second).second;
0260     map_hits[key] = std::pair<double, double>((itr->second).first, energy);
0261 #ifdef EDM_ML_DEBUG
0262     ++i;
0263 #endif
0264   }
0265
0266 #ifdef EDM_ML_DEBUG
0267   edm::LogVerbatim("HGCSim") << "analyzeSimHits: Finds " << map_hits.size() << " hits "
0268                              << " from the Hit Vector of size " << hits.size() << " for type " << type;
0269 #endif
0270   for (const auto& itr : map_hits) {
0271     uint32_t id = (itr.first).first;
0272     double time = (itr.second).first;
0273     double energy = (itr.second).second;
0274     if (type == 0) {
0275       simHitCellIdEE_.push_back(id);
0276       simHitCellEnEE_.push_back(energy);
0277       simHitCellTmEE_.push_back(time);
0278     } else {
0279       simHitCellIdBeam_.push_back(id);
0280       simHitCellEnBeam_.push_back(energy);
0281       simHitCellTmBeam_.push_back(time);
0282     }
0283 #ifdef EDM_ML_DEBUG
0284     edm::LogVerbatim("HGCSim") << "SimHit::ID: " << std::hex << id << std::dec << " T: " << time << " E: " << energy;
0285 #endif
0286   }
0287 }
0288
0289 void HGCalTimingAnalyzer::analyzeSimTracks(edm::Handle<edm::SimTrackContainer> const& SimTk,
0290                                            edm::Handle<edm::SimVertexContainer> const& SimVtx) {
0291   xBeam_ = yBeam_ = zBeam_ = pBeam_ = -1000000;
0292   int vertIndex(-1);
0293   for (edm::SimTrackContainer::const_iterator simTrkItr = SimTk->begin(); simTrkItr != SimTk->end(); simTrkItr++) {
0294 #ifdef EDM_ML_DEBUG
0295     edm::LogVerbatim("HGCSim") << "Track " << simTrkItr->trackId() << " Vertex " << simTrkItr->vertIndex() << " Type "
0296                                << simTrkItr->type() << " Charge " << simTrkItr->charge() << " momentum "
0297                                << simTrkItr->momentum() << " " << simTrkItr->momentum().P();
0298 #endif
0299     if (vertIndex == -1) {
0300       vertIndex = simTrkItr->vertIndex();
0301       pBeam_ = simTrkItr->momentum().P();
0302     }
0303   }
0304   if (vertIndex != -1 && vertIndex < (int)SimVtx->size()) {
0305     edm::SimVertexContainer::const_iterator simVtxItr = SimVtx->begin();
0306     for (int iv = 0; iv < vertIndex; iv++)
0307       simVtxItr++;
0308 #ifdef EDM_ML_DEBUG
0309     edm::LogVerbatim("HGCSim") << "Vertex " << vertIndex << " position " << simVtxItr->position();
0310 #endif
0311     xBeam_ = simVtxItr->position().X();
0312     yBeam_ = simVtxItr->position().Y();
0313     zBeam_ = simVtxItr->position().Z();
0314   }
0315 }
0316
0317 // define this as a plug-in
0318
0319 DEFINE_FWK_MODULE(HGCalTimingAnalyzer);