Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:20:45 

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   // retrieve trigger cells
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   // retrieve clusters
0153   edm::Handle<l1t::HGCalMulticlusterBxCollection> multiclusters_h;
0154   e.getByToken(multiclusters_token_, multiclusters_h);
0155   const l1t::HGCalMulticlusterBxCollection& multiclusters = *multiclusters_h;
0156 
0157   // sim hit association
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   // Associate cells to clusters
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     // loop on 2D clusters inside 3D clusters
0179     for (const auto& cl_ptr : mcl_itr->constituents()) {
0180       // loop on TC inside 2D clusters
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       // Filter cells not included in a multicluster, if requested
0199       if (filter_cells_in_multiclusters_ && mcl_id == 0)
0200         continue;
0201       tc_n_++;
0202       // hardware data
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       // physical values
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       // Links between TC and clusters
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         // search for simhit for all the cells inside the trigger cell
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   // EE
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   //  FH
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   //  BH
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 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!