Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-10-08 22:26:34

0001 #include <memory>
0002 
0003 #include "FWCore/Framework/interface/Frameworkfwd.h"
0004 #include "FWCore/Framework/interface/stream/EDProducer.h"
0005 #include "FWCore/Framework/interface/Event.h"
0006 #include "FWCore/Framework/interface/MakerMacros.h"
0007 #include "FWCore/Framework/interface/ESHandle.h"
0008 #include "FWCore/Framework/interface/EventSetup.h"
0009 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0010 #include "FWCore/Utilities/interface/StreamID.h"
0011 #include "DataFormats/Common/interface/AssociationMap.h"
0012 #include "DataFormats/Common/interface/OneToMany.h"
0013 #include "DataFormats/Common/interface/Ref.h"
0014 #include "DataFormats/L1THGCal/interface/HGCalTriggerCell.h"
0015 #include "DataFormats/L1THGCal/interface/HGCalMulticluster.h"
0016 #include "DataFormats/ForwardDetId/interface/HGCalDetId.h"
0017 #include "SimDataFormats/CaloAnalysis/interface/CaloParticleFwd.h"
0018 #include "SimDataFormats/CaloAnalysis/interface/CaloParticle.h"
0019 #include "SimDataFormats/CaloAnalysis/interface/SimCluster.h"
0020 #include "L1Trigger/L1THGCal/interface/HGCalTriggerGeometryBase.h"
0021 #include "L1Trigger/L1THGCal/interface/backend/HGCalShowerShape.h"
0022 #include "L1Trigger/L1THGCal/interface/backend/HGCalClusteringDummyImpl.h"
0023 #include "L1Trigger/L1THGCal/interface/HGCalTriggerTools.h"
0024 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0025 
0026 class CaloTruthCellsProducer : public edm::stream::EDProducer<> {
0027 public:
0028   explicit CaloTruthCellsProducer(edm::ParameterSet const&);
0029   ~CaloTruthCellsProducer() override;
0030 
0031   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0032 
0033 private:
0034   void beginRun(const edm::Run&, const edm::EventSetup&) override;
0035   void produce(edm::Event&, edm::EventSetup const&) override;
0036 
0037   std::unordered_map<uint32_t, double> makeHitMap(edm::Event const&,
0038                                                   edm::EventSetup const&,
0039                                                   HGCalTriggerGeometryBase const&) const;
0040 
0041   typedef edm::AssociationMap<edm::OneToMany<CaloParticleCollection, l1t::HGCalTriggerCellBxCollection>> CaloToCellsMap;
0042 
0043   bool makeCellsCollection_;
0044   edm::EDGetTokenT<CaloParticleCollection> caloParticlesToken_;
0045   edm::EDGetTokenT<l1t::HGCalTriggerCellBxCollection> triggerCellsToken_;
0046   edm::EDGetTokenT<std::vector<PCaloHit>> simHitsTokenEE_;
0047   edm::EDGetTokenT<std::vector<PCaloHit>> simHitsTokenHEfront_;
0048   edm::EDGetTokenT<std::vector<PCaloHit>> simHitsTokenHEback_;
0049   edm::ESGetToken<HGCalTriggerGeometryBase, CaloGeometryRecord> triggerGeomToken_;
0050   edm::ESHandle<HGCalTriggerGeometryBase> triggerGeomHandle_;
0051 
0052   HGCalClusteringDummyImpl dummyClustering_;
0053   HGCalShowerShape showerShape_;
0054 
0055   HGCalTriggerTools triggerTools_;
0056 };
0057 
0058 CaloTruthCellsProducer::CaloTruthCellsProducer(edm::ParameterSet const& config)
0059     : makeCellsCollection_(config.getParameter<bool>("makeCellsCollection")),
0060       caloParticlesToken_(consumes<CaloParticleCollection>(config.getParameter<edm::InputTag>("caloParticles"))),
0061       triggerCellsToken_(
0062           consumes<l1t::HGCalTriggerCellBxCollection>(config.getParameter<edm::InputTag>("triggerCells"))),
0063       simHitsTokenEE_(consumes<std::vector<PCaloHit>>(config.getParameter<edm::InputTag>("simHitsEE"))),
0064       simHitsTokenHEfront_(consumes<std::vector<PCaloHit>>(config.getParameter<edm::InputTag>("simHitsHEfront"))),
0065       simHitsTokenHEback_(consumes<std::vector<PCaloHit>>(config.getParameter<edm::InputTag>("simHitsHEback"))),
0066       triggerGeomToken_(esConsumes<HGCalTriggerGeometryBase, CaloGeometryRecord, edm::Transition::BeginRun>()),
0067       dummyClustering_(config.getParameterSet("dummyClustering")) {
0068   produces<CaloToCellsMap>();
0069   produces<l1t::HGCalClusterBxCollection>();
0070   produces<l1t::HGCalMulticlusterBxCollection>();
0071   if (makeCellsCollection_)
0072     produces<l1t::HGCalTriggerCellBxCollection>();
0073 }
0074 
0075 CaloTruthCellsProducer::~CaloTruthCellsProducer() {}
0076 
0077 void CaloTruthCellsProducer::beginRun(const edm::Run& /*run*/, const edm::EventSetup& es) {
0078   triggerGeomHandle_ = es.getHandle(triggerGeomToken_);
0079 }
0080 
0081 void CaloTruthCellsProducer::produce(edm::Event& event, edm::EventSetup const& setup) {
0082   edm::Handle<CaloParticleCollection> caloParticlesHandle;
0083   event.getByToken(caloParticlesToken_, caloParticlesHandle);
0084   auto const& caloParticles(*caloParticlesHandle);
0085 
0086   edm::Handle<l1t::HGCalTriggerCellBxCollection> triggerCellsHandle;
0087   event.getByToken(triggerCellsToken_, triggerCellsHandle);
0088   auto const& triggerCells(*triggerCellsHandle);
0089 
0090   auto const& geometry(*triggerGeomHandle_);
0091 
0092   dummyClustering_.setGeometry(triggerGeomHandle_.product());
0093   showerShape_.setGeometry(triggerGeomHandle_.product());
0094   triggerTools_.setGeometry(triggerGeomHandle_.product());
0095 
0096   std::unordered_map<uint32_t, CaloParticleRef> tcToCalo;
0097 
0098   std::unordered_map<uint32_t, double> hitToEnergy(makeHitMap(event, setup, geometry));  // cellId -> sim energy
0099   std::unordered_map<uint32_t, std::pair<double, double>>
0100       tcToEnergies;  // tcId -> {total sim energy, fractioned sim energy}
0101 
0102   for (auto const& he : hitToEnergy) {
0103     DetId hitId(he.first);
0104     // this line will throw if (for whatever reason) hitId is not mapped to a trigger cell id
0105     uint32_t tcId(geometry.getTriggerCellFromCell(hitId));
0106     tcToEnergies[tcId].first += he.second;
0107   }
0108 
0109   // used later to order multiclusters
0110   std::map<int, CaloParticleRef> orderedCaloRefs;
0111 
0112   for (unsigned iP(0); iP != caloParticles.size(); ++iP) {
0113     auto const& caloParticle(caloParticles.at(iP));
0114     if (caloParticle.g4Tracks().at(0).eventId().event() != 0)  // pileup
0115       continue;
0116 
0117     CaloParticleRef ref(caloParticlesHandle, iP);
0118 
0119     SimClusterRefVector const& simClusters(caloParticle.simClusters());
0120     for (auto const& simCluster : simClusters) {
0121       for (auto const& hAndF : simCluster->hits_and_fractions()) {
0122         DetId hitId(hAndF.first);
0123         uint32_t tcId;
0124         try {
0125           tcId = geometry.getTriggerCellFromCell(hitId);
0126         } catch (cms::Exception const& ex) {
0127           edm::LogError("CaloTruthCellsProducer") << ex.what();
0128           continue;
0129         }
0130 
0131         tcToCalo.emplace(tcId, ref);
0132         tcToEnergies[tcId].second += hitToEnergy[hAndF.first] * hAndF.second;
0133       }
0134     }
0135 
0136     // ordered by the gen particle index
0137     int genIndex(caloParticle.g4Tracks().at(0).genpartIndex() - 1);
0138     orderedCaloRefs[genIndex] = ref;
0139   }
0140 
0141   auto outMap(std::make_unique<CaloToCellsMap>(caloParticlesHandle, triggerCellsHandle));
0142   std::unique_ptr<l1t::HGCalTriggerCellBxCollection> outCollection;
0143   if (makeCellsCollection_)
0144     outCollection = std::make_unique<l1t::HGCalTriggerCellBxCollection>();
0145 
0146   typedef edm::Ptr<l1t::HGCalTriggerCell> TriggerCellPtr;
0147   typedef edm::Ptr<l1t::HGCalCluster> ClusterPtr;
0148 
0149   // ClusteringDummyImpl only considers BX 0, so we dump all cells to one vector
0150   std::vector<TriggerCellPtr> triggerCellPtrs;
0151 
0152   // loop through all bunch crossings
0153   for (int bx(triggerCells.getFirstBX()); bx <= triggerCells.getLastBX(); ++bx) {
0154     for (auto cItr(triggerCells.begin(bx)); cItr != triggerCells.end(bx); ++cItr) {
0155       auto const& cell(*cItr);
0156 
0157       auto mapElem(tcToCalo.find(cell.detId()));
0158       if (mapElem == tcToCalo.end())
0159         continue;
0160 
0161       outMap->insert(mapElem->second,
0162                      edm::Ref<l1t::HGCalTriggerCellBxCollection>(triggerCellsHandle, triggerCells.key(cItr)));
0163 
0164       if (makeCellsCollection_) {
0165         auto const& simEnergies(tcToEnergies.at(cell.detId()));
0166         if (simEnergies.first > 0.) {
0167           double eFraction(simEnergies.second / simEnergies.first);
0168 
0169           outCollection->push_back(bx, cell);
0170           auto& newCell((*outCollection)[outCollection->size() - 1]);
0171 
0172           newCell.setMipPt(cell.mipPt() * eFraction);
0173           newCell.setP4(cell.p4() * eFraction);
0174         }
0175       }
0176 
0177       triggerCellPtrs.emplace_back(triggerCellsHandle, triggerCells.key(cItr));
0178     }
0179   }
0180 
0181   event.put(std::move(outMap));
0182   if (makeCellsCollection_)
0183     event.put(std::move(outCollection));
0184 
0185   auto outClusters(std::make_unique<l1t::HGCalClusterBxCollection>());
0186 
0187   auto sortCellPtrs(
0188       [](TriggerCellPtr const& lhs, TriggerCellPtr const& rhs) -> bool { return lhs->mipPt() > rhs->mipPt(); });
0189 
0190   std::sort(triggerCellPtrs.begin(), triggerCellPtrs.end(), sortCellPtrs);
0191   dummyClustering_.clusterizeDummy(triggerCellPtrs, *outClusters);
0192 
0193   std::unordered_map<unsigned, std::vector<unsigned>> caloToClusterIndices;
0194   for (unsigned iC(0); iC != outClusters->size(); ++iC) {
0195     auto const& cluster((*outClusters)[iC]);
0196     // cluster detId and cell detId are identical
0197     auto caloRef(tcToCalo.at(cluster.detId()));
0198     caloToClusterIndices[caloRef.key()].push_back(iC);
0199   }
0200 
0201   auto clustersHandle(event.put(std::move(outClusters)));
0202 
0203   auto outMulticlusters(std::make_unique<l1t::HGCalMulticlusterBxCollection>());
0204 
0205   for (auto const& ocr : orderedCaloRefs) {
0206     auto const& ref(ocr.second);
0207 
0208     if (ref.isNull())  // shouldn't happen
0209       continue;
0210 
0211     auto const& caloParticle(*ref);
0212 
0213     l1t::HGCalMulticluster multicluster;
0214 
0215     for (unsigned iC : caloToClusterIndices[ref.key()]) {
0216       ClusterPtr clPtr(clustersHandle, iC);
0217       multicluster.addConstituent(clPtr, true, 1.);
0218     }
0219 
0220     // Set the gen particle index as the DetId
0221     multicluster.setDetId(caloParticle.g4Tracks().at(0).genpartIndex() - 1);
0222 
0223     auto const& centre(multicluster.centre());
0224     math::PtEtaPhiMLorentzVector multiclusterP4(multicluster.sumPt(), centre.eta(), centre.phi(), 0.);
0225     multicluster.setP4(multiclusterP4);
0226 
0227     showerShape_.fillShapes(multicluster, geometry);
0228 
0229     // not setting the quality flag
0230     // multicluster.setHwQual(id_->decision(multicluster));
0231     // fill H/E
0232     multicluster.saveHOverE();
0233 
0234     outMulticlusters->push_back(0, multicluster);
0235   }
0236 
0237   event.put(std::move(outMulticlusters));
0238 }
0239 
0240 std::unordered_map<uint32_t, double> CaloTruthCellsProducer::makeHitMap(
0241     edm::Event const& event, edm::EventSetup const& setup, HGCalTriggerGeometryBase const& geometry) const {
0242   std::unordered_map<uint32_t, double> hitMap;  // cellId -> sim energy
0243 
0244   typedef std::function<DetId(DetId const&)> DetIdMapper;
0245   typedef std::pair<edm::EDGetTokenT<std::vector<PCaloHit>> const*, DetIdMapper> SimHitSpec;
0246 
0247   SimHitSpec specs[3] = {{&simHitsTokenEE_,
0248                           [this, &geometry](DetId const& simId) -> DetId {
0249                             return this->triggerTools_.simToReco(simId, geometry.eeTopology());
0250                           }},
0251                          {&simHitsTokenHEfront_,
0252                           [this, &geometry](DetId const& simId) -> DetId {
0253                             return this->triggerTools_.simToReco(simId, geometry.fhTopology());
0254                           }},
0255                          {&simHitsTokenHEback_, [this, &geometry](DetId const& simId) -> DetId {
0256                             return this->triggerTools_.simToReco(simId, geometry.hscTopology());
0257                           }}};
0258 
0259   for (auto const& tt : specs) {
0260     edm::Handle<std::vector<PCaloHit>> handle;
0261     event.getByToken(*tt.first, handle);
0262     auto const& simhits(*handle);
0263 
0264     for (auto const& simhit : simhits)
0265       hitMap.emplace(tt.second(simhit.id()), simhit.energy());
0266   }
0267 
0268   return hitMap;
0269 }
0270 
0271 void CaloTruthCellsProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0272   //The following says we do not know what parameters are allowed so do no validation
0273   // Please change this to state exactly what you do use, even if it is no parameters
0274   edm::ParameterSetDescription desc;
0275   desc.setUnknown();
0276   descriptions.addDefault(desc);
0277 }
0278 
0279 DEFINE_FWK_MODULE(CaloTruthCellsProducer);