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
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
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
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
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
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 }
