File indexing completed on 2024-09-07 04:36:56
0001
0002 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0003 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0004 #include "SimDataFormats/CaloTest/interface/HGCalTestNumbering.h"
0005 #include "Geometry/HcalCommonData/interface/HcalHitRelabeller.h"
0006 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0007 #include "DataFormats/L1THGCal/interface/HGCalTriggerCell.h"
0008 #include "DataFormats/L1THGCal/interface/HGCalMulticluster.h"
0009 #include "DataFormats/Common/interface/AssociationMap.h"
0010 #include "DataFormats/Common/interface/OneToMany.h"
0011 #include "SimDataFormats/CaloAnalysis/interface/CaloParticleFwd.h"
0012 #include "SimDataFormats/CaloAnalysis/interface/CaloParticle.h"
0013 #include "DataFormats/ForwardDetId/interface/HGCalTriggerDetId.h"
0014 #include "L1Trigger/L1THGCal/interface/HGCalTriggerGeometryBase.h"
0015 #include "L1Trigger/L1THGCalUtilities/interface/HGCalTriggerNtupleBase.h"
0016 #include "L1Trigger/L1THGCal/interface/HGCalTriggerTools.h"
0017
0018 class HGCalTriggerNtupleHGCTriggerCells : public HGCalTriggerNtupleBase {
0019 public:
0020 HGCalTriggerNtupleHGCTriggerCells(const edm::ParameterSet& conf);
0021 ~HGCalTriggerNtupleHGCTriggerCells() override {}
0022 void initialize(TTree&, const edm::ParameterSet&, edm::ConsumesCollector&&) final;
0023 void fill(const edm::Event& e, const HGCalTriggerNtupleEventSetup& es) final;
0024
0025 private:
0026 double calibrate(double, int, unsigned);
0027 void simhits(const edm::Event& e,
0028 std::unordered_map<uint32_t, double>& simhits_ee,
0029 std::unordered_map<uint32_t, double>& simhits_fh,
0030 std::unordered_map<uint32_t, double>& simhits_bh);
0031 void clear() final;
0032
0033 HGCalTriggerTools triggerTools_;
0034
0035 edm::EDGetToken trigger_cells_token_, multiclusters_token_;
0036 edm::EDGetToken simhits_ee_token_, simhits_fh_token_, simhits_bh_token_;
0037 edm::EDGetToken caloparticles_map_token_;
0038 bool fill_simenergy_;
0039 bool fill_truthmap_;
0040 bool filter_cells_in_multiclusters_;
0041 double keV2fC_;
0042 std::vector<double> fcPerMip_;
0043 std::vector<double> layerWeights_;
0044 std::vector<double> thicknessCorrections_;
0045
0046 int tc_n_;
0047 std::vector<uint32_t> tc_id_;
0048 std::vector<int> tc_subdet_;
0049 std::vector<int> tc_side_;
0050 std::vector<int> tc_layer_;
0051 std::vector<int> tc_waferu_;
0052 std::vector<int> tc_waferv_;
0053 std::vector<int> tc_wafertype_;
0054 std::vector<int> tc_cellu_;
0055 std::vector<int> tc_cellv_;
0056 std::vector<uint32_t> tc_data_;
0057 std::vector<uint32_t> tc_uncompressedCharge_;
0058 std::vector<uint32_t> tc_compressedCharge_;
0059 std::vector<float> tc_mipPt_;
0060 std::vector<float> tc_pt_;
0061 std::vector<float> tc_energy_;
0062 std::vector<float> tc_simenergy_;
0063 std::vector<float> tc_eta_;
0064 std::vector<float> tc_phi_;
0065 std::vector<float> tc_x_;
0066 std::vector<float> tc_y_;
0067 std::vector<float> tc_z_;
0068 std::vector<uint32_t> tc_cluster_id_;
0069 std::vector<uint32_t> tc_multicluster_id_;
0070 std::vector<float> tc_multicluster_pt_;
0071 std::vector<int> tc_genparticle_index_;
0072
0073 typedef edm::AssociationMap<edm::OneToMany<CaloParticleCollection, l1t::HGCalTriggerCellBxCollection>> CaloToCellsMap;
0074 };
0075
0076 DEFINE_EDM_PLUGIN(HGCalTriggerNtupleFactory, HGCalTriggerNtupleHGCTriggerCells, "HGCalTriggerNtupleHGCTriggerCells");
0077
0078 HGCalTriggerNtupleHGCTriggerCells::HGCalTriggerNtupleHGCTriggerCells(const edm::ParameterSet& conf)
0079 : HGCalTriggerNtupleBase(conf),
0080 fill_simenergy_(conf.getParameter<bool>("FillSimEnergy")),
0081 fill_truthmap_(conf.getParameter<bool>("FillTruthMap")),
0082 filter_cells_in_multiclusters_(conf.getParameter<bool>("FilterCellsInMulticlusters")),
0083 keV2fC_(conf.getParameter<double>("keV2fC")),
0084 fcPerMip_(conf.getParameter<std::vector<double>>("fcPerMip")),
0085 layerWeights_(conf.getParameter<std::vector<double>>("layerWeights")),
0086 thicknessCorrections_(conf.getParameter<std::vector<double>>("thicknessCorrections")) {
0087 accessEventSetup_ = false;
0088 }
0089
0090 void HGCalTriggerNtupleHGCTriggerCells::initialize(TTree& tree,
0091 const edm::ParameterSet& conf,
0092 edm::ConsumesCollector&& collector) {
0093 trigger_cells_token_ =
0094 collector.consumes<l1t::HGCalTriggerCellBxCollection>(conf.getParameter<edm::InputTag>("TriggerCells"));
0095 multiclusters_token_ =
0096 collector.consumes<l1t::HGCalMulticlusterBxCollection>(conf.getParameter<edm::InputTag>("Multiclusters"));
0097
0098 if (fill_simenergy_) {
0099 simhits_ee_token_ = collector.consumes<edm::PCaloHitContainer>(conf.getParameter<edm::InputTag>("eeSimHits"));
0100 simhits_fh_token_ = collector.consumes<edm::PCaloHitContainer>(conf.getParameter<edm::InputTag>("fhSimHits"));
0101 simhits_bh_token_ = collector.consumes<edm::PCaloHitContainer>(conf.getParameter<edm::InputTag>("bhSimHits"));
0102 }
0103
0104 if (fill_truthmap_)
0105 caloparticles_map_token_ =
0106 collector.consumes<CaloToCellsMap>(conf.getParameter<edm::InputTag>("caloParticlesToCells"));
0107
0108 std::string prefix(conf.getUntrackedParameter<std::string>("Prefix", "tc"));
0109
0110 std::string bname;
0111 auto withPrefix([&prefix, &bname](char const* vname) -> char const* {
0112 bname = prefix + "_" + vname;
0113 return bname.c_str();
0114 });
0115
0116 tree.Branch(withPrefix("n"), &tc_n_, (prefix + "_n/I").c_str());
0117 tree.Branch(withPrefix("id"), &tc_id_);
0118 tree.Branch(withPrefix("subdet"), &tc_subdet_);
0119 tree.Branch(withPrefix("zside"), &tc_side_);
0120 tree.Branch(withPrefix("layer"), &tc_layer_);
0121 tree.Branch(withPrefix("waferu"), &tc_waferu_);
0122 tree.Branch(withPrefix("waferv"), &tc_waferv_);
0123 tree.Branch(withPrefix("wafertype"), &tc_wafertype_);
0124 tree.Branch(withPrefix("cellu"), &tc_cellu_);
0125 tree.Branch(withPrefix("cellv"), &tc_cellv_);
0126 tree.Branch(withPrefix("data"), &tc_data_);
0127 tree.Branch(withPrefix("uncompressedCharge"), &tc_uncompressedCharge_);
0128 tree.Branch(withPrefix("compressedCharge"), &tc_compressedCharge_);
0129 tree.Branch(withPrefix("pt"), &tc_pt_);
0130 tree.Branch(withPrefix("mipPt"), &tc_mipPt_);
0131 tree.Branch(withPrefix("energy"), &tc_energy_);
0132 if (fill_simenergy_)
0133 tree.Branch(withPrefix("simenergy"), &tc_simenergy_);
0134 tree.Branch(withPrefix("eta"), &tc_eta_);
0135 tree.Branch(withPrefix("phi"), &tc_phi_);
0136 tree.Branch(withPrefix("x"), &tc_x_);
0137 tree.Branch(withPrefix("y"), &tc_y_);
0138 tree.Branch(withPrefix("z"), &tc_z_);
0139 tree.Branch(withPrefix("cluster_id"), &tc_cluster_id_);
0140 tree.Branch(withPrefix("multicluster_id"), &tc_multicluster_id_);
0141 tree.Branch(withPrefix("multicluster_pt"), &tc_multicluster_pt_);
0142 if (fill_truthmap_)
0143 tree.Branch(withPrefix("genparticle_index"), &tc_genparticle_index_);
0144 }
0145
0146 void HGCalTriggerNtupleHGCTriggerCells::fill(const edm::Event& e, const HGCalTriggerNtupleEventSetup& es) {
0147
0148 edm::Handle<l1t::HGCalTriggerCellBxCollection> trigger_cells_h;
0149 e.getByToken(trigger_cells_token_, trigger_cells_h);
0150 const l1t::HGCalTriggerCellBxCollection& trigger_cells = *trigger_cells_h;
0151
0152
0153 edm::Handle<l1t::HGCalMulticlusterBxCollection> multiclusters_h;
0154 e.getByToken(multiclusters_token_, multiclusters_h);
0155 const l1t::HGCalMulticlusterBxCollection& multiclusters = *multiclusters_h;
0156
0157
0158 std::unordered_map<uint32_t, double> simhits_ee;
0159 std::unordered_map<uint32_t, double> simhits_fh;
0160 std::unordered_map<uint32_t, double> simhits_bh;
0161 if (fill_simenergy_)
0162 simhits(e, simhits_ee, simhits_fh, simhits_bh);
0163
0164 edm::Handle<CaloToCellsMap> caloparticles_map_h;
0165 std::unordered_map<uint32_t, unsigned> cell_to_genparticle;
0166 if (fill_truthmap_) {
0167 e.getByToken(caloparticles_map_token_, caloparticles_map_h);
0168 for (auto& keyval : *caloparticles_map_h) {
0169 for (auto& tcref : keyval.val)
0170 cell_to_genparticle.emplace(tcref->detId(), keyval.key->g4Tracks().at(0).genpartIndex() - 1);
0171 }
0172 }
0173
0174
0175 std::unordered_map<uint32_t, uint32_t> cell2cluster;
0176 std::unordered_map<uint32_t, l1t::HGCalMulticlusterBxCollection::const_iterator> cell2multicluster;
0177 for (auto mcl_itr = multiclusters.begin(0); mcl_itr != multiclusters.end(0); mcl_itr++) {
0178
0179 for (const auto& cl_ptr : mcl_itr->constituents()) {
0180
0181 for (const auto& tc_ptr : cl_ptr.second->constituents()) {
0182 cell2cluster.emplace(tc_ptr.second->detId(), cl_ptr.second->detId());
0183 cell2multicluster.emplace(tc_ptr.second->detId(), mcl_itr);
0184 }
0185 }
0186 }
0187
0188 triggerTools_.setGeometry(es.geometry.product());
0189
0190 clear();
0191 for (auto tc_itr = trigger_cells.begin(0); tc_itr != trigger_cells.end(0); tc_itr++) {
0192 if (tc_itr->hwPt() > 0) {
0193 auto cl_itr = cell2cluster.find(tc_itr->detId());
0194 auto mcl_itr = cell2multicluster.find(tc_itr->detId());
0195 uint32_t cl_id = (cl_itr != cell2cluster.end() ? cl_itr->second : 0);
0196 uint32_t mcl_id = (mcl_itr != cell2multicluster.end() ? mcl_itr->second->detId() : 0);
0197 float mcl_pt = (mcl_itr != cell2multicluster.end() ? mcl_itr->second->pt() : 0.);
0198
0199 if (filter_cells_in_multiclusters_ && mcl_id == 0)
0200 continue;
0201 tc_n_++;
0202
0203 DetId id(tc_itr->detId());
0204 tc_id_.emplace_back(tc_itr->detId());
0205 tc_side_.emplace_back(triggerTools_.zside(id));
0206 tc_layer_.emplace_back(triggerTools_.layerWithOffset(id));
0207 if (id.det() == DetId::HGCalTrigger) {
0208 HGCalTriggerDetId idtrg(id);
0209 tc_subdet_.emplace_back(idtrg.subdet());
0210 tc_waferu_.emplace_back(idtrg.waferU());
0211 tc_waferv_.emplace_back(idtrg.waferV());
0212 tc_wafertype_.emplace_back(idtrg.type());
0213 tc_cellu_.emplace_back(idtrg.triggerCellU());
0214 tc_cellv_.emplace_back(idtrg.triggerCellV());
0215 } else if (id.det() == DetId::HGCalHSc) {
0216 HGCScintillatorDetId idsci(id);
0217 tc_subdet_.emplace_back(idsci.subdet());
0218 tc_waferu_.emplace_back(-999);
0219 tc_waferv_.emplace_back(-999);
0220 tc_wafertype_.emplace_back(idsci.type());
0221 tc_cellu_.emplace_back(idsci.ietaAbs());
0222 tc_cellv_.emplace_back(idsci.iphi());
0223 } else {
0224 throw cms::Exception("InvalidHGCalTriggerDetid")
0225 << "Found unexpected trigger cell detid to be filled in HGCal Trigger Cell ntuple.";
0226 }
0227 tc_data_.emplace_back(tc_itr->hwPt());
0228 tc_uncompressedCharge_.emplace_back(tc_itr->uncompressedCharge());
0229 tc_compressedCharge_.emplace_back(tc_itr->compressedCharge());
0230 tc_mipPt_.emplace_back(tc_itr->mipPt());
0231
0232 tc_pt_.emplace_back(tc_itr->pt());
0233 tc_energy_.emplace_back(tc_itr->energy());
0234 tc_eta_.emplace_back(tc_itr->eta());
0235 tc_phi_.emplace_back(tc_itr->phi());
0236 tc_x_.emplace_back(tc_itr->position().x());
0237 tc_y_.emplace_back(tc_itr->position().y());
0238 tc_z_.emplace_back(tc_itr->position().z());
0239
0240 tc_cluster_id_.emplace_back(cl_id);
0241 tc_multicluster_id_.emplace_back(mcl_id);
0242 tc_multicluster_pt_.emplace_back(mcl_pt);
0243
0244 if (fill_simenergy_) {
0245 double energy = 0;
0246 unsigned layer = triggerTools_.layerWithOffset(id);
0247
0248 for (uint32_t c_id : triggerTools_.getTriggerGeometry()->getCellsFromTriggerCell(id)) {
0249 int thickness = triggerTools_.thicknessIndex(c_id);
0250 if (triggerTools_.isEm(id)) {
0251 auto itr = simhits_ee.find(c_id);
0252 if (itr != simhits_ee.end())
0253 energy += calibrate(itr->second, thickness, layer);
0254 } else if (triggerTools_.isSilicon(id)) {
0255 auto itr = simhits_fh.find(c_id);
0256 if (itr != simhits_fh.end())
0257 energy += calibrate(itr->second, thickness, layer);
0258 } else {
0259 auto itr = simhits_bh.find(c_id);
0260 if (itr != simhits_bh.end())
0261 energy += itr->second;
0262 }
0263 }
0264 tc_simenergy_.emplace_back(energy);
0265 }
0266 }
0267
0268 if (fill_truthmap_) {
0269 auto itr(cell_to_genparticle.find(tc_itr->detId()));
0270 if (itr == cell_to_genparticle.end())
0271 tc_genparticle_index_.push_back(-1);
0272 else
0273 tc_genparticle_index_.push_back(itr->second);
0274 }
0275 }
0276 }
0277
0278 double HGCalTriggerNtupleHGCTriggerCells::calibrate(double energy, int thickness, unsigned layer) {
0279 double fcPerMip = fcPerMip_[thickness];
0280 double thicknessCorrection = thicknessCorrections_[thickness];
0281 double layerWeight = layerWeights_[layer];
0282 double TeV2GeV = 1.e3;
0283 return energy * keV2fC_ / fcPerMip * layerWeight * TeV2GeV / thicknessCorrection;
0284 }
0285
0286 void HGCalTriggerNtupleHGCTriggerCells::simhits(const edm::Event& e,
0287 std::unordered_map<uint32_t, double>& simhits_ee,
0288 std::unordered_map<uint32_t, double>& simhits_fh,
0289 std::unordered_map<uint32_t, double>& simhits_bh) {
0290 edm::Handle<edm::PCaloHitContainer> ee_simhits_h;
0291 e.getByToken(simhits_ee_token_, ee_simhits_h);
0292 const edm::PCaloHitContainer& ee_simhits = *ee_simhits_h;
0293 edm::Handle<edm::PCaloHitContainer> fh_simhits_h;
0294 e.getByToken(simhits_fh_token_, fh_simhits_h);
0295 const edm::PCaloHitContainer& fh_simhits = *fh_simhits_h;
0296 edm::Handle<edm::PCaloHitContainer> bh_simhits_h;
0297 e.getByToken(simhits_bh_token_, bh_simhits_h);
0298 const edm::PCaloHitContainer& bh_simhits = *bh_simhits_h;
0299
0300
0301 for (const auto& simhit : ee_simhits) {
0302 DetId id = triggerTools_.simToReco(simhit.id(), triggerTools_.getTriggerGeometry()->eeTopology());
0303 if (id.rawId() == 0)
0304 continue;
0305 auto itr_insert = simhits_ee.emplace(id, 0.);
0306 itr_insert.first->second += simhit.energy();
0307 }
0308
0309 for (const auto& simhit : fh_simhits) {
0310 DetId id = triggerTools_.simToReco(simhit.id(), triggerTools_.getTriggerGeometry()->fhTopology());
0311 if (id.rawId() == 0)
0312 continue;
0313 auto itr_insert = simhits_fh.emplace(id, 0.);
0314 itr_insert.first->second += simhit.energy();
0315 }
0316
0317 for (const auto& simhit : bh_simhits) {
0318 DetId id = triggerTools_.simToReco(simhit.id(), triggerTools_.getTriggerGeometry()->hscTopology());
0319 if (id.rawId() == 0)
0320 continue;
0321 auto itr_insert = simhits_bh.emplace(id, 0.);
0322 itr_insert.first->second += simhit.energy();
0323 }
0324 }
0325
0326 void HGCalTriggerNtupleHGCTriggerCells::clear() {
0327 tc_n_ = 0;
0328 tc_id_.clear();
0329 tc_subdet_.clear();
0330 tc_side_.clear();
0331 tc_layer_.clear();
0332 tc_waferu_.clear();
0333 tc_waferv_.clear();
0334 tc_wafertype_.clear();
0335 tc_cellu_.clear();
0336 tc_cellv_.clear();
0337 tc_data_.clear();
0338 tc_uncompressedCharge_.clear();
0339 tc_compressedCharge_.clear();
0340 tc_mipPt_.clear();
0341 tc_pt_.clear();
0342 tc_energy_.clear();
0343 tc_simenergy_.clear();
0344 tc_eta_.clear();
0345 tc_phi_.clear();
0346 tc_x_.clear();
0347 tc_y_.clear();
0348 tc_z_.clear();
0349 tc_cluster_id_.clear();
0350 tc_multicluster_id_.clear();
0351 tc_multicluster_pt_.clear();
0352 tc_genparticle_index_.clear();
0353 }