L1THGCalUtilities/plugins/CaloTruthCellsProducer.cc

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