Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-04-27 00:16:57

0001 #ifndef __RecoLocalCalo_HGCRecProducers_HGCalLayerClusterProducer_H__
0002 #define __RecoLocalCalo_HGCRecProducers_HGCalLayerClusterProducer_H__
0003 
0004 // user include files
0005 #include "FWCore/Framework/interface/Frameworkfwd.h"
0006 #include "FWCore/Framework/interface/stream/EDProducer.h"
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008 
0009 #include "FWCore/Framework/interface/Event.h"
0010 #include "FWCore/Framework/interface/ESHandle.h"
0011 #include "FWCore/Framework/interface/MakerMacros.h"
0012 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0013 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0014 #include "FWCore/ParameterSet/interface/PluginDescription.h"
0015 
0016 #include "RecoParticleFlow/PFClusterProducer/interface/RecHitTopologicalCleanerBase.h"
0017 #include "RecoParticleFlow/PFClusterProducer/interface/SeedFinderBase.h"
0018 #include "RecoParticleFlow/PFClusterProducer/interface/InitialClusteringStepBase.h"
0019 #include "RecoParticleFlow/PFClusterProducer/interface/PFClusterBuilderBase.h"
0020 #include "RecoParticleFlow/PFClusterProducer/interface/PFCPositionCalculatorBase.h"
0021 #include "RecoParticleFlow/PFClusterProducer/interface/PFClusterEnergyCorrectorBase.h"
0022 #include "RecoLocalCalo/HGCalRecProducers/interface/ComputeClusterTime.h"
0023 
0024 #include "RecoLocalCalo/HGCalRecProducers/interface/HGCalLayerClusterAlgoFactory.h"
0025 #include "RecoLocalCalo/HGCalRecAlgos/interface/HGCalDepthPreClusterer.h"
0026 
0027 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0028 #include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
0029 
0030 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
0031 #include "DataFormats/Common/interface/ValueMap.h"
0032 
0033 using Density = hgcal_clustering::Density;
0034 
0035 class HGCalLayerClusterProducer : public edm::stream::EDProducer<> {
0036 public:
0037   HGCalLayerClusterProducer(const edm::ParameterSet&);
0038   ~HGCalLayerClusterProducer() override {}
0039   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0040 
0041   void produce(edm::Event&, const edm::EventSetup&) override;
0042 
0043 private:
0044   edm::EDGetTokenT<HGCRecHitCollection> hits_ee_token;
0045   edm::EDGetTokenT<HGCRecHitCollection> hits_fh_token;
0046   edm::EDGetTokenT<HGCRecHitCollection> hits_bh_token;
0047   edm::EDGetTokenT<HGCRecHitCollection> hits_hfnose_token;
0048 
0049   reco::CaloCluster::AlgoId algoId;
0050 
0051   std::unique_ptr<HGCalClusteringAlgoBase> algo;
0052   bool doSharing;
0053   std::string detector;
0054 
0055   std::string timeClname;
0056   double timeOffset;
0057   unsigned int nHitsTime;
0058 };
0059 
0060 DEFINE_FWK_MODULE(HGCalLayerClusterProducer);
0061 
0062 HGCalLayerClusterProducer::HGCalLayerClusterProducer(const edm::ParameterSet& ps)
0063     : algoId(reco::CaloCluster::undefined),
0064       doSharing(ps.getParameter<bool>("doSharing")),
0065       detector(ps.getParameter<std::string>("detector")),  // one of EE, FH, BH or "all"
0066       timeClname(ps.getParameter<std::string>("timeClname")),
0067       timeOffset(ps.getParameter<double>("timeOffset")),
0068       nHitsTime(ps.getParameter<unsigned int>("nHitsTime")) {
0069   if (detector == "HFNose") {
0070     hits_hfnose_token = consumes<HGCRecHitCollection>(ps.getParameter<edm::InputTag>("HFNoseInput"));
0071     algoId = reco::CaloCluster::hfnose;
0072   } else if (detector == "all") {
0073     hits_ee_token = consumes<HGCRecHitCollection>(ps.getParameter<edm::InputTag>("HGCEEInput"));
0074     hits_fh_token = consumes<HGCRecHitCollection>(ps.getParameter<edm::InputTag>("HGCFHInput"));
0075     hits_bh_token = consumes<HGCRecHitCollection>(ps.getParameter<edm::InputTag>("HGCBHInput"));
0076     algoId = reco::CaloCluster::hgcal_mixed;
0077   } else if (detector == "EE") {
0078     hits_ee_token = consumes<HGCRecHitCollection>(ps.getParameter<edm::InputTag>("HGCEEInput"));
0079     algoId = reco::CaloCluster::hgcal_em;
0080   } else if (detector == "FH") {
0081     hits_fh_token = consumes<HGCRecHitCollection>(ps.getParameter<edm::InputTag>("HGCFHInput"));
0082     algoId = reco::CaloCluster::hgcal_had;
0083   } else {
0084     hits_bh_token = consumes<HGCRecHitCollection>(ps.getParameter<edm::InputTag>("HGCBHInput"));
0085     algoId = reco::CaloCluster::hgcal_had;
0086   }
0087 
0088   auto pluginPSet = ps.getParameter<edm::ParameterSet>("plugin");
0089   if (detector == "HFNose") {
0090     algo = HGCalLayerClusterAlgoFactory::get()->create("HFNoseCLUE", pluginPSet, consumesCollector());
0091     algo->setAlgoId(algoId, true);
0092   } else {
0093     algo = HGCalLayerClusterAlgoFactory::get()->create(
0094         pluginPSet.getParameter<std::string>("type"), pluginPSet, consumesCollector());
0095     algo->setAlgoId(algoId);
0096   }
0097 
0098   produces<std::vector<float>>("InitialLayerClustersMask");
0099   produces<std::vector<reco::BasicCluster>>();
0100   produces<std::vector<reco::BasicCluster>>("sharing");
0101   //density
0102   produces<Density>();
0103   //time for layer clusters
0104   produces<edm::ValueMap<std::pair<float, float>>>(timeClname);
0105 }
0106 
0107 void HGCalLayerClusterProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0108   // hgcalLayerClusters
0109   edm::ParameterSetDescription desc;
0110   edm::ParameterSetDescription pluginDesc;
0111   pluginDesc.addNode(edm::PluginDescription<HGCalLayerClusterAlgoFactory>("type", "CLUE", true));
0112 
0113   desc.add<edm::ParameterSetDescription>("plugin", pluginDesc);
0114   desc.add<std::string>("detector", "all")
0115       ->setComment("all (does not include HFNose); other options: EE, FH, HFNose; other value defaults to BH");
0116   desc.add<bool>("doSharing", false);
0117   desc.add<edm::InputTag>("HFNoseInput", edm::InputTag("HGCalRecHit", "HGCHFNoseRecHits"));
0118   desc.add<edm::InputTag>("HGCEEInput", edm::InputTag("HGCalRecHit", "HGCEERecHits"));
0119   desc.add<edm::InputTag>("HGCFHInput", edm::InputTag("HGCalRecHit", "HGCHEFRecHits"));
0120   desc.add<edm::InputTag>("HGCBHInput", edm::InputTag("HGCalRecHit", "HGCHEBRecHits"));
0121   desc.add<std::string>("timeClname", "timeLayerCluster");
0122   desc.add<double>("timeOffset", 0.0);
0123   desc.add<unsigned int>("nHitsTime", 3);
0124   descriptions.add("hgcalLayerClusters", desc);
0125 }
0126 
0127 void HGCalLayerClusterProducer::produce(edm::Event& evt, const edm::EventSetup& es) {
0128   edm::Handle<HGCRecHitCollection> hfnose_hits;
0129   edm::Handle<HGCRecHitCollection> ee_hits;
0130   edm::Handle<HGCRecHitCollection> fh_hits;
0131   edm::Handle<HGCRecHitCollection> bh_hits;
0132 
0133   std::unique_ptr<std::vector<reco::BasicCluster>> clusters(new std::vector<reco::BasicCluster>),
0134       clusters_sharing(new std::vector<reco::BasicCluster>);
0135   auto density = std::make_unique<Density>();
0136 
0137   algo->getEventSetup(es);
0138 
0139   //make a map detid-rechit
0140   // NB for the moment just host EE and FH hits
0141   // timing in digi for BH not implemented for now
0142   std::unordered_map<uint32_t, const HGCRecHit*> hitmap;
0143 
0144   switch (algoId) {
0145     case reco::CaloCluster::hfnose:
0146       evt.getByToken(hits_hfnose_token, hfnose_hits);
0147       algo->populate(*hfnose_hits);
0148       for (auto const& it : *hfnose_hits)
0149         hitmap[it.detid().rawId()] = &(it);
0150       break;
0151     case reco::CaloCluster::hgcal_em:
0152       evt.getByToken(hits_ee_token, ee_hits);
0153       algo->populate(*ee_hits);
0154       for (auto const& it : *ee_hits)
0155         hitmap[it.detid().rawId()] = &(it);
0156       break;
0157     case reco::CaloCluster::hgcal_had:
0158       evt.getByToken(hits_fh_token, fh_hits);
0159       evt.getByToken(hits_bh_token, bh_hits);
0160       if (fh_hits.isValid()) {
0161         algo->populate(*fh_hits);
0162         for (auto const& it : *fh_hits)
0163           hitmap[it.detid().rawId()] = &(it);
0164       } else if (bh_hits.isValid()) {
0165         algo->populate(*bh_hits);
0166       }
0167       break;
0168     case reco::CaloCluster::hgcal_mixed:
0169       evt.getByToken(hits_ee_token, ee_hits);
0170       algo->populate(*ee_hits);
0171       for (auto const& it : *ee_hits) {
0172         hitmap[it.detid().rawId()] = &(it);
0173       }
0174       evt.getByToken(hits_fh_token, fh_hits);
0175       algo->populate(*fh_hits);
0176       for (auto const& it : *fh_hits) {
0177         hitmap[it.detid().rawId()] = &(it);
0178       }
0179       evt.getByToken(hits_bh_token, bh_hits);
0180       algo->populate(*bh_hits);
0181       break;
0182     default:
0183       break;
0184   }
0185   algo->makeClusters();
0186   *clusters = algo->getClusters(false);
0187   if (doSharing)
0188     *clusters_sharing = algo->getClusters(true);
0189 
0190   auto clusterHandle = evt.put(std::move(clusters));
0191   auto clusterHandleSharing = evt.put(std::move(clusters_sharing), "sharing");
0192 
0193   //Keep the density
0194   *density = algo->getDensity();
0195   evt.put(std::move(density));
0196 
0197   edm::PtrVector<reco::BasicCluster> clusterPtrs, clusterPtrsSharing;
0198 
0199   std::vector<std::pair<float, float>> times;
0200   times.reserve(clusterHandle->size());
0201 
0202   for (unsigned i = 0; i < clusterHandle->size(); ++i) {
0203     edm::Ptr<reco::BasicCluster> ptr(clusterHandle, i);
0204     clusterPtrs.push_back(ptr);
0205 
0206     std::pair<float, float> timeCl(-99., -1.);
0207 
0208     const reco::CaloCluster& sCl = (*clusterHandle)[i];
0209     if (sCl.size() >= nHitsTime) {
0210       std::vector<float> timeClhits;
0211       std::vector<float> timeErrorClhits;
0212 
0213       for (auto const& hit : sCl.hitsAndFractions()) {
0214         auto finder = hitmap.find(hit.first);
0215         if (finder == hitmap.end())
0216           continue;
0217 
0218         //time is computed wrt  0-25ns + offset and set to -1 if no time
0219         const HGCRecHit* rechit = finder->second;
0220         float rhTimeE = rechit->timeError();
0221         //check on timeError to exclude scintillator
0222         if (rhTimeE < 0.)
0223           continue;
0224         timeClhits.push_back(rechit->time() - timeOffset);
0225         timeErrorClhits.push_back(1. / (rhTimeE * rhTimeE));
0226       }
0227       hgcalsimclustertime::ComputeClusterTime timeEstimator;
0228       timeCl = timeEstimator.fixSizeHighestDensity(timeClhits, timeErrorClhits, nHitsTime);
0229     }
0230     times.push_back(timeCl);
0231   }
0232   std::unique_ptr<std::vector<float>> layerClustersMask(new std::vector<float>);
0233   layerClustersMask->resize(clusterHandle->size(), 1.0);
0234   evt.put(std::move(layerClustersMask), "InitialLayerClustersMask");
0235 
0236   auto timeCl = std::make_unique<edm::ValueMap<std::pair<float, float>>>();
0237   edm::ValueMap<std::pair<float, float>>::Filler filler(*timeCl);
0238   filler.insert(clusterHandle, times.begin(), times.end());
0239   filler.fill();
0240   evt.put(std::move(timeCl), timeClname);
0241 
0242   if (doSharing) {
0243     for (unsigned i = 0; i < clusterHandleSharing->size(); ++i) {
0244       edm::Ptr<reco::BasicCluster> ptr(clusterHandleSharing, i);
0245       clusterPtrsSharing.push_back(ptr);
0246     }
0247   }
0248   algo->reset();
0249 }
0250 
0251 #endif  //__RecoLocalCalo_HGCRecProducers_HGCalLayerClusterProducer_H__