PFClusterProducer/plugins/BarrelLayerClusterProducer.cc

0001 #ifndef __RecoParticleFlow_PFClusterProducer_BarrelLayerClusterProducer_H__
0002 #define __RecoParticleFlow_PFClusterProducer_BarrelLayerClusterProducer_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 "RecoParticleFlow/PFClusterProducer/interface/CaloRecHitResolutionProvider.h"
0023
0024 #include "RecoLocalCalo/HGCalRecProducers/interface/ComputeClusterTime.h"
0025
0026 #include "RecoLocalCalo/HGCalRecProducers/interface/HGCalLayerClusterAlgoFactory.h"
0027 #include "RecoLocalCalo/HGCalRecAlgos/interface/HGCalDepthPreClusterer.h"
0028
0029 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0030 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0031
0032 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
0033 #include "DataFormats/ParticleFlowReco/interface/PFRecHit.h"
0034
0035 #include "DataFormats/Common/interface/ValueMap.h"
0036
0037 #include "RecoParticleFlow/PFClusterProducer/plugins/BarrelCLUEAlgo.h"
0038 #include "RecoLocalCalo/HGCalRecProducers/interface/BarrelTilesConstants.h"
0039
0040 using Density = hgcal_clustering::Density;
0041
0042 class BarrelLayerClusterProducer : public edm::stream::EDProducer<> {
0043 public:
0044   BarrelLayerClusterProducer(const edm::ParameterSet&);
0045   ~BarrelLayerClusterProducer() override {}
0046   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0047
0048   void produce(edm::Event&, const edm::EventSetup&) override;
0049
0050 private:
0051   reco::CaloCluster::AlgoId algoId_;
0052   std::string timeClname_;
0053   unsigned int nHitsTime_;
0054   edm::EDGetTokenT<reco::PFRecHitCollection> hits_token_;
0055   edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeomToken_;
0056
0057   std::unique_ptr<HGCalClusteringAlgoBase> algo_;
0058   hgcal::RecHitTools rhtools_;
0059   std::unique_ptr<CaloRecHitResolutionProvider> timeResolutionCalc_;
0060
0061   void setAlgoId(std::string& type);
0062 };
0063
0064 DEFINE_FWK_MODULE(BarrelLayerClusterProducer);
0065
0066 BarrelLayerClusterProducer::BarrelLayerClusterProducer(const edm::ParameterSet& ps)
0067     : algoId_(reco::CaloCluster::undefined),
0068       timeClname_(ps.getParameter<std::string>("timeClname")),
0069       nHitsTime_(ps.getParameter<unsigned int>("nHitsTime")),
0070       hits_token_{consumes<reco::PFRecHitCollection>(ps.getParameter<edm::InputTag>("recHits"))},
0071       caloGeomToken_{consumesCollector().esConsumes<CaloGeometry, CaloGeometryRecord>()} {
0072   auto pluginPSet = ps.getParameter<edm::ParameterSet>("plugin");
0073   std::string type = pluginPSet.getParameter<std::string>("type");
0074   algo_ = HGCalLayerClusterAlgoFactory::get()->create(type, pluginPSet);
0075   setAlgoId(type);
0076
0077   algo_->setThresholds(consumesCollector().esConsumes<EcalPFRecHitThresholds, EcalPFRecHitThresholdsRcd>(),
0078                        consumesCollector().esConsumes<HcalPFCuts, HcalPFCutsRcd>());
0079
0080   timeResolutionCalc_ = std::make_unique<CaloRecHitResolutionProvider>(ps.getParameterSet("timeResolutionCalc"));
0081   produces<std::vector<float>>("InitialLayerClustersMask");
0082   produces<std::vector<reco::BasicCluster>>();
0083   produces<edm::ValueMap<std::pair<float, float>>>(timeClname_);
0084 }
0085
0086 void BarrelLayerClusterProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0087   edm::ParameterSetDescription desc;
0088
0089   edm::ParameterSetDescription pluginDesc;
0090   pluginDesc.addNode(edm::PluginDescription<HGCalLayerClusterAlgoFactory>("type", "EBCLUE", true));
0091
0092   desc.add<edm::ParameterSetDescription>("plugin", pluginDesc);
0093
0094   desc.add<edm::InputTag>("recHits", edm::InputTag("particleFlowRecHitECAL", ""));
0095
0096   edm::ParameterSetDescription timeResolutionCalcDesc;
0097   timeResolutionCalcDesc.addNode(edm::ParameterDescription<double>("noiseTerm", 1.10889, true) and
0098                                  edm::ParameterDescription<double>("constantTerm", 0.428192, true) and
0099                                  edm::ParameterDescription<double>("corrTermLowE", 0.0510871, true) and
0100                                  edm::ParameterDescription<double>("threshLowE", 0.5, true) and
0101                                  edm::ParameterDescription<double>("constantTermLowE", 0.0, true) and
0102                                  edm::ParameterDescription<double>("noiseTermLowE", 1.31883, true) and
0103                                  edm::ParameterDescription<double>("threshHighE", 5.0, true));
0104   desc.add<edm::ParameterSetDescription>("timeResolutionCalc", timeResolutionCalcDesc);
0105
0106   desc.add<std::string>("timeClname", "timeLayerCluster");
0107   desc.add<unsigned int>("nHitsTime", 3);
0108   descriptions.add("barrelLayerClusters", desc);
0109 }
0110
0111 void BarrelLayerClusterProducer::produce(edm::Event& evt, const edm::EventSetup& es) {
0112   edm::Handle<reco::PFRecHitCollection> hits;
0113   edm::ESHandle<CaloGeometry> geom = es.getHandle(caloGeomToken_);
0114   rhtools_.setGeometry(*geom);
0115
0116   //auto density = std::make_unique<Density>();
0117   algo_->getEventSetup(es, rhtools_);
0118
0119   //make a map detid-rechit
0120   // NB for the moment just host EE and FH hits
0121   // timing in digi for BH not implemented for now
0122   std::unordered_map<uint32_t, const reco::PFRecHit*> hitmap;
0123
0124   evt.getByToken(hits_token_, hits);
0125   algo_->populate(*hits);
0126   for (auto& hit : *hits) {
0127     hitmap[hit.detId()] = &(hit);
0128   }
0129   algo_->makeClusters();
0130
0131   std::unique_ptr<std::vector<reco::CaloCluster>> clusters(new std::vector<reco::CaloCluster>);
0132   *clusters = algo_->getClusters(false);
0133   auto clusterHandle = evt.put(std::move(clusters));
0134
0135   edm::PtrVector<reco::BasicCluster> clusterPtrs;  //, clusterPtrsSharing;
0136
0137   std::vector<std::pair<float, float>> times;
0138   times.reserve(clusterHandle->size());
0139
0140   for (unsigned i = 0; i < clusterHandle->size(); ++i) {
0141     edm::Ptr<reco::BasicCluster> ptr(clusterHandle, i);
0142     clusterPtrs.push_back(ptr);
0143
0144     std::pair<float, float> timeCl(-99., -1.);
0145
0146     const reco::CaloCluster& sCl = (*clusterHandle)[i];
0147     if (sCl.size() >= nHitsTime_) {
0148       std::vector<float> timeClhits;
0149       std::vector<float> timeErrorClhits;
0150
0151       for (auto const& hit : sCl.hitsAndFractions()) {
0152         auto finder = hitmap.find(hit.first);
0153         if (finder == hitmap.end())
0154           continue;
0155
0156         //time is computed wrt  0-25ns + offset and set to -1 if no time
0157         const reco::PFRecHit* rechit = finder->second;
0158         float rhTimeE = timeResolutionCalc_->timeResolution2(rechit->energy());
0159         if (rhTimeE < 0.)
0160           continue;
0161         timeClhits.push_back(rechit->time());
0162         timeErrorClhits.push_back(1. / rhTimeE);
0163       }
0164       hgcalsimclustertime::ComputeClusterTime timeEstimator;
0165       timeCl = timeEstimator.fixSizeHighestDensity(timeClhits, timeErrorClhits, nHitsTime_);
0166     }
0167     times.push_back(timeCl);
0168   }
0169
0170   std::unique_ptr<std::vector<float>> layerClustersMask(new std::vector<float>);
0171   layerClustersMask->resize(clusterHandle->size(), 1.0);
0172   evt.put(std::move(layerClustersMask), "InitialLayerClustersMask");
0173
0174   auto timeCl = std::make_unique<edm::ValueMap<std::pair<float, float>>>();
0175   edm::ValueMap<std::pair<float, float>>::Filler filler(*timeCl);
0176   filler.insert(clusterHandle, times.begin(), times.end());
0177   filler.fill();
0178   evt.put(std::move(timeCl), timeClname_);
0179
0180   algo_->reset();
0181 }
0182
0183 void BarrelLayerClusterProducer::setAlgoId(std::string& type) {
0184   if (type == "EBCLUE") {
0185     algoId_ = reco::CaloCluster::barrel_em;
0186   } else if (type == "HBCLUE") {
0187     algoId_ = reco::CaloCluster::barrel_had;
0188   } else {
0189     throw cms::Exception("InvalidPlugin") << "Invalid plugin type: " << type << std::endl;
0190   }
0191 }
0192
0193 #endif  //__RecoLocalCalo_HGCRecProducers_BarrelLayerClusterProducer_H__