Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-10-25 09:55:18

0001 #include "DataFormats/HGCDigi/interface/HGCDigiCollections.h"
0002 #include "L1Trigger/L1THGCalUtilities/interface/HGCalTriggerNtupleBase.h"
0003 #include "L1Trigger/L1THGCal/interface/HGCalTriggerGeometryBase.h"
0004 #include "FWCore/Framework/interface/ESHandle.h"
0005 
0006 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0007 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0008 
0009 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0010 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0011 #include "L1Trigger/L1THGCal/interface/HGCalTriggerTools.h"
0012 
0013 class HGCalTriggerNtupleHGCDigis : public HGCalTriggerNtupleBase {
0014 public:
0015   HGCalTriggerNtupleHGCDigis(const edm::ParameterSet& conf);
0016   ~HGCalTriggerNtupleHGCDigis() override{};
0017   void initialize(TTree&, const edm::ParameterSet&, edm::ConsumesCollector&&) final;
0018   void fill(const edm::Event& e, const HGCalTriggerNtupleEventSetup& es) final;
0019 
0020 private:
0021   void simhits(const edm::Event& e,
0022                std::unordered_map<uint32_t, double>& simhits_ee,
0023                std::unordered_map<uint32_t, double>& simhits_fh,
0024                std::unordered_map<uint32_t, double>& simhits_bh);
0025   void clear() final;
0026 
0027   edm::EDGetToken ee_token_, fh_token_, bh_token_;
0028   bool is_Simhit_comp_;
0029   edm::EDGetToken SimHits_inputee_, SimHits_inputfh_, SimHits_inputbh_;
0030 
0031   std::vector<unsigned int> digiBXselect_;
0032   static constexpr unsigned kDigiSize_ = 5;
0033 
0034   HGCalTriggerTools triggerTools_;
0035 
0036   int hgcdigi_n_;
0037   std::vector<int> hgcdigi_id_;
0038   std::vector<int> hgcdigi_subdet_;
0039   std::vector<int> hgcdigi_side_;
0040   std::vector<int> hgcdigi_layer_;
0041   std::vector<int> hgcdigi_wafertype_;
0042   std::vector<float> hgcdigi_eta_;
0043   std::vector<float> hgcdigi_phi_;
0044   std::vector<float> hgcdigi_z_;
0045   std::vector<std::vector<uint32_t>> hgcdigi_data_;
0046   std::vector<std::vector<int>> hgcdigi_isadc_;
0047   std::vector<float> hgcdigi_simenergy_;
0048   std::vector<int> hgcdigi_waferu_;
0049   std::vector<int> hgcdigi_waferv_;
0050   std::vector<int> hgcdigi_cellu_;
0051   std::vector<int> hgcdigi_cellv_;
0052 
0053   int bhdigi_n_;
0054   std::vector<int> bhdigi_id_;
0055   std::vector<int> bhdigi_subdet_;
0056   std::vector<int> bhdigi_side_;
0057   std::vector<int> bhdigi_layer_;
0058   std::vector<int> bhdigi_ieta_;
0059   std::vector<int> bhdigi_iphi_;
0060   std::vector<float> bhdigi_eta_;
0061   std::vector<float> bhdigi_phi_;
0062   std::vector<float> bhdigi_z_;
0063   std::vector<std::vector<uint32_t>> bhdigi_data_;
0064   std::vector<std::vector<int>> bhdigi_isadc_;
0065   std::vector<float> bhdigi_simenergy_;
0066 };
0067 
0068 DEFINE_EDM_PLUGIN(HGCalTriggerNtupleFactory, HGCalTriggerNtupleHGCDigis, "HGCalTriggerNtupleHGCDigis");
0069 
0070 HGCalTriggerNtupleHGCDigis::HGCalTriggerNtupleHGCDigis(const edm::ParameterSet& conf) : HGCalTriggerNtupleBase(conf) {
0071   accessEventSetup_ = false;
0072   is_Simhit_comp_ = conf.getParameter<bool>("isSimhitComp");
0073   digiBXselect_ = conf.getParameter<std::vector<unsigned int>>("digiBXselect");
0074 
0075   if (digiBXselect_.empty()) {
0076     throw cms::Exception("BadInitialization") << "digiBXselect vector is empty";
0077   }
0078   if (*std::max_element(digiBXselect_.begin(), digiBXselect_.end()) >= kDigiSize_) {
0079     throw cms::Exception("BadInitialization")
0080         << "digiBXselect vector requests a BX outside of maximum size of digis (" << kDigiSize_ << " BX)";
0081   }
0082   //sort and check for duplicates
0083   std::sort(digiBXselect_.begin(), digiBXselect_.end());
0084   if (std::unique(digiBXselect_.begin(), digiBXselect_.end()) != digiBXselect_.end()) {
0085     throw cms::Exception("BadInitialization") << "digiBXselect vector contains duplicate BX values";
0086   }
0087 }
0088 
0089 void HGCalTriggerNtupleHGCDigis::initialize(TTree& tree,
0090                                             const edm::ParameterSet& conf,
0091                                             edm::ConsumesCollector&& collector) {
0092   ee_token_ = collector.consumes<HGCalDigiCollection>(conf.getParameter<edm::InputTag>("HGCDigisEE"));
0093   fh_token_ = collector.consumes<HGCalDigiCollection>(conf.getParameter<edm::InputTag>("HGCDigisFH"));
0094   bh_token_ = collector.consumes<HGCalDigiCollection>(conf.getParameter<edm::InputTag>("HGCDigisBH"));
0095   if (is_Simhit_comp_) {
0096     SimHits_inputee_ = collector.consumes<edm::PCaloHitContainer>(conf.getParameter<edm::InputTag>("eeSimHits"));
0097     SimHits_inputfh_ = collector.consumes<edm::PCaloHitContainer>(conf.getParameter<edm::InputTag>("fhSimHits"));
0098     SimHits_inputbh_ = collector.consumes<edm::PCaloHitContainer>(conf.getParameter<edm::InputTag>("bhSimHits"));
0099   }
0100 
0101   hgcdigi_data_.resize(digiBXselect_.size());
0102   hgcdigi_isadc_.resize(digiBXselect_.size());
0103   bhdigi_data_.resize(digiBXselect_.size());
0104   bhdigi_isadc_.resize(digiBXselect_.size());
0105 
0106   tree.Branch("hgcdigi_n", &hgcdigi_n_, "hgcdigi_n/I");
0107   tree.Branch("hgcdigi_id", &hgcdigi_id_);
0108   tree.Branch("hgcdigi_subdet", &hgcdigi_subdet_);
0109   tree.Branch("hgcdigi_zside", &hgcdigi_side_);
0110   tree.Branch("hgcdigi_layer", &hgcdigi_layer_);
0111   tree.Branch("hgcdigi_wafertype", &hgcdigi_wafertype_);
0112   tree.Branch("hgcdigi_eta", &hgcdigi_eta_);
0113   tree.Branch("hgcdigi_phi", &hgcdigi_phi_);
0114   tree.Branch("hgcdigi_z", &hgcdigi_z_);
0115   std::string bname;
0116   auto withBX([&bname](char const* vname, unsigned int bx) -> char const* {
0117     bname = std::string(vname) + "_BX" + to_string(bx);
0118     return bname.c_str();
0119   });
0120   for (unsigned int i = 0; i < digiBXselect_.size(); i++) {
0121     unsigned int bxi = digiBXselect_[i];
0122     tree.Branch(withBX("hgcdigi_data", bxi), &hgcdigi_data_[i]);
0123     tree.Branch(withBX("hgcdigi_isadc", bxi), &hgcdigi_isadc_[i]);
0124   }
0125   tree.Branch("hgcdigi_waferu", &hgcdigi_waferu_);
0126   tree.Branch("hgcdigi_waferv", &hgcdigi_waferv_);
0127   tree.Branch("hgcdigi_cellu", &hgcdigi_cellu_);
0128   tree.Branch("hgcdigi_cellv", &hgcdigi_cellv_);
0129   if (is_Simhit_comp_)
0130     tree.Branch("hgcdigi_simenergy", &hgcdigi_simenergy_);
0131 
0132   tree.Branch("bhdigi_n", &bhdigi_n_, "bhdigi_n/I");
0133   tree.Branch("bhdigi_id", &bhdigi_id_);
0134   tree.Branch("bhdigi_subdet", &bhdigi_subdet_);
0135   tree.Branch("bhdigi_zside", &bhdigi_side_);
0136   tree.Branch("bhdigi_layer", &bhdigi_layer_);
0137   tree.Branch("bhdigi_ieta", &bhdigi_ieta_);
0138   tree.Branch("bhdigi_iphi", &bhdigi_iphi_);
0139   tree.Branch("bhdigi_eta", &bhdigi_eta_);
0140   tree.Branch("bhdigi_phi", &bhdigi_phi_);
0141   tree.Branch("bhdigi_z", &bhdigi_z_);
0142   for (unsigned int i = 0; i < digiBXselect_.size(); i++) {
0143     unsigned int bxi = digiBXselect_[i];
0144     tree.Branch(withBX("bhdigi_data", bxi), &bhdigi_data_[i]);
0145     tree.Branch(withBX("bhdigi_isadc", bxi), &bhdigi_isadc_[i]);
0146   }
0147   if (is_Simhit_comp_)
0148     tree.Branch("bhdigi_simenergy", &bhdigi_simenergy_);
0149 }
0150 
0151 void HGCalTriggerNtupleHGCDigis::fill(const edm::Event& e, const HGCalTriggerNtupleEventSetup& es) {
0152   edm::Handle<HGCalDigiCollection> ee_digis_h;
0153   e.getByToken(ee_token_, ee_digis_h);
0154   const HGCalDigiCollection& ee_digis = *ee_digis_h;
0155   edm::Handle<HGCalDigiCollection> fh_digis_h;
0156   e.getByToken(fh_token_, fh_digis_h);
0157   const HGCalDigiCollection& fh_digis = *fh_digis_h;
0158   edm::Handle<HGCalDigiCollection> bh_digis_h;
0159   e.getByToken(bh_token_, bh_digis_h);
0160   const HGCalDigiCollection& bh_digis = *bh_digis_h;
0161 
0162   triggerTools_.setGeometry(es.geometry.product());
0163 
0164   // sim hit association
0165   std::unordered_map<uint32_t, double> simhits_ee;
0166   std::unordered_map<uint32_t, double> simhits_fh;
0167   std::unordered_map<uint32_t, double> simhits_bh;
0168   if (is_Simhit_comp_)
0169     simhits(e, simhits_ee, simhits_fh, simhits_bh);
0170 
0171   clear();
0172   hgcdigi_n_ = ee_digis.size() + fh_digis.size();
0173   hgcdigi_id_.reserve(hgcdigi_n_);
0174   hgcdigi_subdet_.reserve(hgcdigi_n_);
0175   hgcdigi_side_.reserve(hgcdigi_n_);
0176   hgcdigi_layer_.reserve(hgcdigi_n_);
0177   hgcdigi_wafertype_.reserve(hgcdigi_n_);
0178   hgcdigi_eta_.reserve(hgcdigi_n_);
0179   hgcdigi_phi_.reserve(hgcdigi_n_);
0180   hgcdigi_z_.reserve(hgcdigi_n_);
0181   for (unsigned int i = 0; i < digiBXselect_.size(); i++) {
0182     hgcdigi_data_[i].reserve(hgcdigi_n_);
0183     hgcdigi_isadc_[i].reserve(hgcdigi_n_);
0184   }
0185   hgcdigi_waferu_.reserve(hgcdigi_n_);
0186   hgcdigi_waferv_.reserve(hgcdigi_n_);
0187   hgcdigi_cellu_.reserve(hgcdigi_n_);
0188   hgcdigi_cellv_.reserve(hgcdigi_n_);
0189   if (is_Simhit_comp_)
0190     hgcdigi_simenergy_.reserve(hgcdigi_n_);
0191 
0192   bhdigi_n_ = bh_digis.size();
0193   bhdigi_id_.reserve(bhdigi_n_);
0194   bhdigi_subdet_.reserve(bhdigi_n_);
0195   bhdigi_side_.reserve(bhdigi_n_);
0196   bhdigi_layer_.reserve(bhdigi_n_);
0197   bhdigi_ieta_.reserve(bhdigi_n_);
0198   bhdigi_iphi_.reserve(bhdigi_n_);
0199   bhdigi_eta_.reserve(bhdigi_n_);
0200   bhdigi_phi_.reserve(bhdigi_n_);
0201   bhdigi_z_.reserve(bhdigi_n_);
0202   for (unsigned int i = 0; i < digiBXselect_.size(); i++) {
0203     bhdigi_data_[i].reserve(bhdigi_n_);
0204     bhdigi_isadc_[i].reserve(bhdigi_n_);
0205   }
0206   if (is_Simhit_comp_)
0207     bhdigi_simenergy_.reserve(bhdigi_n_);
0208 
0209   for (const auto& digi : ee_digis) {
0210     const DetId id(digi.id());
0211     hgcdigi_id_.emplace_back(id.rawId());
0212     hgcdigi_subdet_.emplace_back(id.det());
0213     hgcdigi_side_.emplace_back(triggerTools_.zside(id));
0214     hgcdigi_layer_.emplace_back(triggerTools_.layerWithOffset(id));
0215     GlobalPoint cellpos = triggerTools_.getTriggerGeometry()->eeGeometry()->getPosition(id.rawId());
0216     hgcdigi_eta_.emplace_back(cellpos.eta());
0217     hgcdigi_phi_.emplace_back(cellpos.phi());
0218     hgcdigi_z_.emplace_back(cellpos.z());
0219     for (unsigned int i = 0; i < digiBXselect_.size(); i++) {
0220       hgcdigi_data_[i].emplace_back(digi[digiBXselect_[i]].data());
0221       hgcdigi_isadc_[i].emplace_back(!digi[digiBXselect_[i]].mode());
0222     }
0223     const HGCSiliconDetId idsi(digi.id());
0224     hgcdigi_waferu_.emplace_back(idsi.waferU());
0225     hgcdigi_waferv_.emplace_back(idsi.waferV());
0226     hgcdigi_wafertype_.emplace_back(idsi.type());
0227     hgcdigi_cellu_.emplace_back(idsi.cellU());
0228     hgcdigi_cellv_.emplace_back(idsi.cellV());
0229     if (is_Simhit_comp_) {
0230       double hit_energy = 0;
0231       auto itr = simhits_ee.find(id);
0232       if (itr != simhits_ee.end())
0233         hit_energy = itr->second;
0234       hgcdigi_simenergy_.emplace_back(hit_energy);
0235     }
0236   }
0237 
0238   for (const auto& digi : fh_digis) {
0239     const DetId id(digi.id());
0240     hgcdigi_id_.emplace_back(id.rawId());
0241     hgcdigi_subdet_.emplace_back(id.det());
0242     hgcdigi_side_.emplace_back(triggerTools_.zside(id));
0243     hgcdigi_layer_.emplace_back(triggerTools_.layerWithOffset(id));
0244     GlobalPoint cellpos = triggerTools_.getTriggerGeometry()->hsiGeometry()->getPosition(id.rawId());
0245     hgcdigi_eta_.emplace_back(cellpos.eta());
0246     hgcdigi_phi_.emplace_back(cellpos.phi());
0247     hgcdigi_z_.emplace_back(cellpos.z());
0248     for (unsigned int i = 0; i < digiBXselect_.size(); i++) {
0249       hgcdigi_data_[i].emplace_back(digi[digiBXselect_[i]].data());
0250       hgcdigi_isadc_[i].emplace_back(!digi[digiBXselect_[i]].mode());
0251     }
0252     const HGCSiliconDetId idsi(digi.id());
0253     hgcdigi_waferu_.emplace_back(idsi.waferU());
0254     hgcdigi_waferv_.emplace_back(idsi.waferV());
0255     hgcdigi_wafertype_.emplace_back(idsi.type());
0256     hgcdigi_cellu_.emplace_back(idsi.cellU());
0257     hgcdigi_cellv_.emplace_back(idsi.cellV());
0258     if (is_Simhit_comp_) {
0259       double hit_energy = 0;
0260       auto itr = simhits_fh.find(id);
0261       if (itr != simhits_fh.end())
0262         hit_energy = itr->second;
0263       hgcdigi_simenergy_.emplace_back(hit_energy);
0264     }
0265   }
0266 
0267   for (const auto& digi : bh_digis) {
0268     const DetId id(digi.id());
0269     bhdigi_id_.emplace_back(id.rawId());
0270     bhdigi_subdet_.emplace_back(id.det());
0271     bhdigi_side_.emplace_back(triggerTools_.zside(id));
0272     bhdigi_layer_.emplace_back(triggerTools_.layerWithOffset(id));
0273     GlobalPoint cellpos = triggerTools_.getTriggerGeometry()->hscGeometry()->getPosition(id.rawId());
0274     bhdigi_eta_.emplace_back(cellpos.eta());
0275     bhdigi_phi_.emplace_back(cellpos.phi());
0276     bhdigi_z_.emplace_back(cellpos.z());
0277     for (unsigned int i = 0; i < digiBXselect_.size(); i++) {
0278       bhdigi_data_[i].emplace_back(digi[digiBXselect_[i]].data());
0279       bhdigi_isadc_[i].emplace_back(!digi[digiBXselect_[i]].mode());
0280     }
0281     const HGCScintillatorDetId idsci(digi.id());
0282     bhdigi_ieta_.emplace_back(idsci.ietaAbs());
0283     bhdigi_iphi_.emplace_back(idsci.iphi());
0284     if (is_Simhit_comp_) {
0285       double hit_energy = 0;
0286       auto itr = simhits_bh.find(id);
0287       if (itr != simhits_bh.end())
0288         hit_energy = itr->second;
0289       bhdigi_simenergy_.emplace_back(hit_energy);
0290     }
0291   }
0292 }
0293 
0294 void HGCalTriggerNtupleHGCDigis::simhits(const edm::Event& e,
0295                                          std::unordered_map<uint32_t, double>& simhits_ee,
0296                                          std::unordered_map<uint32_t, double>& simhits_fh,
0297                                          std::unordered_map<uint32_t, double>& simhits_bh) {
0298   edm::Handle<edm::PCaloHitContainer> ee_simhits_h;
0299   e.getByToken(SimHits_inputee_, ee_simhits_h);
0300   const edm::PCaloHitContainer& ee_simhits = *ee_simhits_h;
0301   edm::Handle<edm::PCaloHitContainer> fh_simhits_h;
0302   e.getByToken(SimHits_inputfh_, fh_simhits_h);
0303   const edm::PCaloHitContainer& fh_simhits = *fh_simhits_h;
0304   edm::Handle<edm::PCaloHitContainer> bh_simhits_h;
0305   e.getByToken(SimHits_inputbh_, bh_simhits_h);
0306   const edm::PCaloHitContainer& bh_simhits = *bh_simhits_h;
0307 
0308   //EE
0309   for (const auto& simhit : ee_simhits) {
0310     DetId id = triggerTools_.simToReco(simhit.id(), triggerTools_.getTriggerGeometry()->eeTopology());
0311     if (id.rawId() == 0)
0312       continue;
0313     auto itr_insert = simhits_ee.emplace(id, 0.);
0314     itr_insert.first->second += simhit.energy();
0315   }
0316   //  FH
0317   for (const auto& simhit : fh_simhits) {
0318     DetId id = triggerTools_.simToReco(simhit.id(), triggerTools_.getTriggerGeometry()->fhTopology());
0319     if (id.rawId() == 0)
0320       continue;
0321     auto itr_insert = simhits_fh.emplace(id, 0.);
0322     itr_insert.first->second += simhit.energy();
0323   }
0324   //  BH
0325   for (const auto& simhit : bh_simhits) {
0326     DetId id = triggerTools_.simToReco(simhit.id(), triggerTools_.getTriggerGeometry()->hscTopology());
0327     if (id.rawId() == 0)
0328       continue;
0329     auto itr_insert = simhits_bh.emplace(id, 0.);
0330     itr_insert.first->second += simhit.energy();
0331   }
0332 }
0333 
0334 void HGCalTriggerNtupleHGCDigis::clear() {
0335   hgcdigi_n_ = 0;
0336   hgcdigi_id_.clear();
0337   hgcdigi_subdet_.clear();
0338   hgcdigi_side_.clear();
0339   hgcdigi_layer_.clear();
0340   hgcdigi_waferu_.clear();
0341   hgcdigi_waferv_.clear();
0342   hgcdigi_wafertype_.clear();
0343   hgcdigi_cellu_.clear();
0344   hgcdigi_cellv_.clear();
0345   hgcdigi_eta_.clear();
0346   hgcdigi_phi_.clear();
0347   hgcdigi_z_.clear();
0348   for (unsigned int i = 0; i < digiBXselect_.size(); i++) {
0349     hgcdigi_data_[i].clear();
0350     hgcdigi_isadc_[i].clear();
0351   }
0352   if (is_Simhit_comp_)
0353     hgcdigi_simenergy_.clear();
0354 
0355   bhdigi_n_ = 0;
0356   bhdigi_id_.clear();
0357   bhdigi_subdet_.clear();
0358   bhdigi_side_.clear();
0359   bhdigi_layer_.clear();
0360   bhdigi_ieta_.clear();
0361   bhdigi_iphi_.clear();
0362   bhdigi_eta_.clear();
0363   bhdigi_phi_.clear();
0364   bhdigi_z_.clear();
0365   for (unsigned int i = 0; i < digiBXselect_.size(); i++) {
0366     bhdigi_data_[i].clear();
0367     bhdigi_isadc_[i].clear();
0368   }
0369   if (is_Simhit_comp_)
0370     bhdigi_simenergy_.clear();
0371 }