HGCalRecProducers/plugins/HGCalLayerClusterProducer.cc

0001 // Authors: Olivie Franklova - olivie.abigail.franklova@cern.ch
0002 // Date: 03/2023
0003 // @file create layer clusters
0004
0005 #define DEBUG_CLUSTERS_ALPAKA 0
0006
0007 // user include files
0008 #include "FWCore/Framework/interface/Frameworkfwd.h"
0009 #include "FWCore/Framework/interface/stream/EDProducer.h"
0010 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0011
0012 #include "FWCore/Framework/interface/Event.h"
0013 #include "FWCore/Framework/interface/ESHandle.h"
0014 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0015 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0016 #include "FWCore/ParameterSet/interface/PluginDescription.h"
0017
0018 #include "RecoParticleFlow/PFClusterProducer/interface/RecHitTopologicalCleanerBase.h"
0019 #include "RecoParticleFlow/PFClusterProducer/interface/SeedFinderBase.h"
0020 #include "RecoParticleFlow/PFClusterProducer/interface/InitialClusteringStepBase.h"
0021 #include "RecoParticleFlow/PFClusterProducer/interface/PFClusterBuilderBase.h"
0022 #include "RecoParticleFlow/PFClusterProducer/interface/PFCPositionCalculatorBase.h"
0023 #include "RecoParticleFlow/PFClusterProducer/interface/PFClusterEnergyCorrectorBase.h"
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/HGCalGeometry/interface/HGCalGeometry.h"
0031
0032 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
0033 #include "DataFormats/Common/interface/ValueMap.h"
0034
0035 #include "FWCore/Framework/interface/ConsumesCollector.h"
0036
0037 #if DEBUG_CLUSTERS_ALPAKA
0038 #include "RecoLocalCalo/HGCalRecProducers/interface/DumpClustersDetails.h"
0039 #endif
0040
0041 class HGCalLayerClusterProducer : public edm::stream::EDProducer<> {
0042 public:
0043   /**
0044    * @brief Constructor with parameter settings - which can be changed in hgcalLayerCluster_cff.py.
0045    * Constructor will set all variables by input param ps.
0046    * algoID variables will be set accordingly to the detector type.
0047    *
0048    * @param[in] ps parametr set to set variables
0049   */
0050   HGCalLayerClusterProducer(const edm::ParameterSet&);
0051   ~HGCalLayerClusterProducer() override {}
0052   /**
0053    * @brief Method fill description which will be used in pyhton file.
0054    *
0055    * @param[out] description to be fill
0056   */
0057   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0058
0059   /**
0060    * @brief Method run the algoritm to get clusters.
0061    *
0062    * @param[in, out] evt from get info and put result
0063    * @param[in] es to get event setup info
0064   */
0065   void produce(edm::Event&, const edm::EventSetup&) override;
0066
0067 private:
0068   edm::EDGetTokenT<HGCRecHitCollection> hits_token_;
0069
0070   reco::CaloCluster::AlgoId algoId_;
0071
0072   std::unique_ptr<HGCalClusteringAlgoBase> algo_;
0073   std::string detector_;
0074
0075   std::string timeClname_;
0076   unsigned int hitsTime_;
0077
0078   // for calculate position
0079   std::vector<double> thresholdW0_;
0080   double positionDeltaRho2_;
0081   hgcal::RecHitTools rhtools_;
0082   edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeomToken_;
0083   const bool calculatePositionInAlgo_;
0084 #if DEBUG_CLUSTERS_ALPAKA
0085   std::string moduleLabel_;
0086 #endif
0087
0088   /**
0089    * @brief Sets algoId accordingly to the detector type
0090   */
0091   void setAlgoId();
0092
0093   /**
0094    * @brief Counts position for all points in the cluster
0095    *
0096    * @param[in] hitmap hitmap to find correct RecHit
0097    * @param[in] hitsAndFraction all hits in the cluster
0098    * @return counted position
0099   */
0100   math::XYZPoint calculatePosition(std::unordered_map<uint32_t, const HGCRecHit*>& hitmap,
0101                                    const std::vector<std::pair<DetId, float>>& hitsAndFractions);
0102
0103   /**
0104    * @brief Counts time for all points in the cluster
0105    *
0106    * @param[in] hitmap hitmap to find correct RecHit only for silicon (not for BH-HSci)
0107    * @param[in] hitsAndFraction all hits in the cluster
0108    * @return counted time
0109   */
0110   std::pair<float, float> calculateTime(std::unordered_map<uint32_t, const HGCRecHit*>& hitmap,
0111                                         const std::vector<std::pair<DetId, float>>& hitsAndFractions,
0112                                         size_t sizeCluster);
0113 };
0114
0115 HGCalLayerClusterProducer::HGCalLayerClusterProducer(const edm::ParameterSet& ps)
0116     : algoId_(reco::CaloCluster::undefined),
0117       detector_(ps.getParameter<std::string>("detector")),  // one of EE, FH, BH, HFNose
0118       timeClname_(ps.getParameter<std::string>("timeClname")),
0119       hitsTime_(ps.getParameter<unsigned int>("nHitsTime")),
0120       caloGeomToken_(consumesCollector().esConsumes<CaloGeometry, CaloGeometryRecord>()),
0121       calculatePositionInAlgo_(ps.getParameter<bool>("calculatePositionInAlgo")) {
0122 #if DEBUG_CLUSTERS_ALPAKA
0123   moduleLabel_ = ps.getParameter<std::string>("@module_label");
0124 #endif
0125   setAlgoId();  //sets algo id according to detector type
0126   hits_token_ = consumes<HGCRecHitCollection>(ps.getParameter<edm::InputTag>("recHits"));
0127
0128   auto pluginPSet = ps.getParameter<edm::ParameterSet>("plugin");
0129   if (detector_ == "HFNose") {
0130     algo_ = HGCalLayerClusterAlgoFactory::get()->create("HFNoseCLUE", pluginPSet);
0131     algo_->setAlgoId(algoId_, true);
0132   } else {
0133     algo_ = HGCalLayerClusterAlgoFactory::get()->create(pluginPSet.getParameter<std::string>("type"), pluginPSet);
0134     algo_->setAlgoId(algoId_);
0135   }
0136   thresholdW0_ = pluginPSet.getParameter<std::vector<double>>("thresholdW0");
0137   positionDeltaRho2_ = pluginPSet.getParameter<double>("positionDeltaRho2");
0138
0139   produces<std::vector<float>>("InitialLayerClustersMask");
0140   produces<std::vector<reco::BasicCluster>>();
0141   //time for layer clusters
0142   produces<edm::ValueMap<std::pair<float, float>>>(timeClname_);
0143 }
0144
0145 void HGCalLayerClusterProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0146   // hgcalLayerClusters
0147   edm::ParameterSetDescription desc;
0148   edm::ParameterSetDescription pluginDesc;
0149   pluginDesc.addNode(edm::PluginDescription<HGCalLayerClusterAlgoFactory>("type", "SiCLUE", true));
0150
0151   desc.add<edm::ParameterSetDescription>("plugin", pluginDesc);
0152   desc.add<std::string>("detector", "EE")->setComment("options EE, FH, BH,  HFNose; other value defaults to EE");
0153   desc.add<edm::InputTag>("recHits", edm::InputTag("HGCalRecHit", "HGCEERecHits"));
0154   desc.add<std::string>("timeClname", "timeLayerCluster");
0155   desc.add<unsigned int>("nHitsTime", 3);
0156   desc.add<bool>("calculatePositionInAlgo", true);
0157   descriptions.add("hgcalLayerClusters", desc);
0158 }
0159
0160 math::XYZPoint HGCalLayerClusterProducer::calculatePosition(
0161     std::unordered_map<uint32_t, const HGCRecHit*>& hitmap,
0162     const std::vector<std::pair<DetId, float>>& hitsAndFractions) {
0163   float total_weight = 0.f;
0164   float maxEnergyValue = 0.f;
0165   DetId maxEnergyIndex;
0166   float x = 0.f;
0167   float y = 0.f;
0168
0169   for (auto const& hit : hitsAndFractions) {
0170     //time is computed wrt  0-25ns + offset and set to -1 if no time
0171     const HGCRecHit* rechit = hitmap[hit.first];
0172     total_weight += rechit->energy();
0173     if (rechit->energy() > maxEnergyValue) {
0174       maxEnergyValue = rechit->energy();
0175       maxEnergyIndex = rechit->detid();
0176     }
0177   }
0178   float total_weight_log = 0.f;
0179   auto thick = rhtools_.getSiThickIndex(maxEnergyIndex);
0180   const GlobalPoint positionMaxEnergy(rhtools_.getPosition(maxEnergyIndex));
0181   for (auto const& hit : hitsAndFractions) {
0182     //time is computed wrt  0-25ns + offset and set to -1 if no time
0183     const HGCRecHit* rechit = hitmap[hit.first];
0184
0185     const GlobalPoint position(rhtools_.getPosition(rechit->detid()));
0186
0187     if (thick != -1) {  //silicon
0188       //for silicon only just use 1+6 cells = 1.3cm for all thicknesses
0189       const float d1 = position.x() - positionMaxEnergy.x();
0190       const float d2 = position.y() - positionMaxEnergy.y();
0191       if ((d1 * d1 + d2 * d2) > positionDeltaRho2_)
0192         continue;
0193
0194       float Wi = std::max(thresholdW0_[thick] + std::log(rechit->energy() / total_weight), 0.);
0195       x += position.x() * Wi;
0196       y += position.y() * Wi;
0197       total_weight_log += Wi;
0198     } else {  //scintillator
0199       x += position.x() * rechit->energy();
0200       y += position.y() * rechit->energy();
0201     }
0202   }
0203   if (thick != -1) {
0204     total_weight = total_weight_log;
0205   }
0206   if (total_weight != 0.) {
0207     float inv_tot_weight = 1.f / total_weight;
0208     return math::XYZPoint(x * inv_tot_weight, y * inv_tot_weight, positionMaxEnergy.z());
0209   } else {
0210     return {positionMaxEnergy.x(), positionMaxEnergy.y(), positionMaxEnergy.z()};
0211   }
0212 }
0213
0214 std::pair<float, float> HGCalLayerClusterProducer::calculateTime(
0215     std::unordered_map<uint32_t, const HGCRecHit*>& hitmap,
0216     const std::vector<std::pair<DetId, float>>& hitsAndFractions,
0217     size_t sizeCluster) {
0218   std::pair<float, float> timeCl(-99., -1.);
0219
0220   if (sizeCluster >= hitsTime_) {
0221     std::vector<float> timeClhits;
0222     std::vector<float> timeErrorClhits;
0223
0224     for (auto const& hit : hitsAndFractions) {
0225       //time is computed wrt  0-25ns + offset and set to -1 if no time
0226       const HGCRecHit* rechit = hitmap[hit.first];
0227
0228       float rhTimeE = rechit->timeError();
0229       //check on timeError to exclude scintillator
0230       if (rhTimeE < 0.f)
0231         continue;
0232       timeClhits.push_back(rechit->time());
0233       timeErrorClhits.push_back(1.f / (rhTimeE * rhTimeE));
0234     }
0235     hgcalsimclustertime::ComputeClusterTime timeEstimator;
0236     timeCl = timeEstimator.fixSizeHighestDensity(timeClhits, timeErrorClhits, hitsTime_);
0237   }
0238   return timeCl;
0239 }
0240 void HGCalLayerClusterProducer::produce(edm::Event& evt, const edm::EventSetup& es) {
0241   edm::Handle<HGCRecHitCollection> hits;
0242
0243   std::unique_ptr<std::vector<reco::BasicCluster>> clusters(new std::vector<reco::BasicCluster>);
0244
0245   edm::ESHandle<CaloGeometry> geom = es.getHandle(caloGeomToken_);
0246   rhtools_.setGeometry(*geom);
0247   algo_->getEventSetup(es, rhtools_);
0248
0249   //make a map detid-rechit
0250   // NB for the moment just host EE and FH hits
0251   // timing in digi for BH not implemented for now
0252   std::unordered_map<uint32_t, const HGCRecHit*> hitmap;
0253
0254   evt.getByToken(hits_token_, hits);
0255   algo_->populate(*hits);
0256   for (auto const& it : *hits) {
0257     hitmap[it.detid().rawId()] = &(it);
0258   }
0259
0260   algo_->makeClusters();
0261   *clusters = algo_->getClusters(false);
0262
0263   std::vector<std::pair<float, float>> times;
0264   times.reserve(clusters->size());
0265
0266   for (unsigned i = 0; i < clusters->size(); ++i) {
0267     reco::CaloCluster& sCl = (*clusters)[i];
0268     if (!calculatePositionInAlgo_) {
0269       sCl.setPosition(calculatePosition(hitmap, sCl.hitsAndFractions()));
0270     }
0271     if (detector_ != "BH") {
0272       times.push_back(calculateTime(hitmap, sCl.hitsAndFractions(), sCl.size()));
0273     } else {
0274       times.push_back(std::pair<float, float>(-99.f, -1.f));
0275     }
0276   }
0277
0278 #if DEBUG_CLUSTERS_ALPAKA
0279   hgcalUtils::DumpClusters dumper;
0280   auto runNumber = evt.eventAuxiliary().run();
0281   auto lumiNumber = evt.eventAuxiliary().luminosityBlock();
0282   auto evtNumber = evt.eventAuxiliary().id().event();
0283
0284   dumper.dumpInfos(*clusters, moduleLabel_, runNumber, lumiNumber, evtNumber, true);
0285 #endif
0286
0287   auto clusterHandle = evt.put(std::move(clusters));
0288
0289   if (detector_ == "HFNose") {
0290     std::unique_ptr<std::vector<float>> layerClustersMask(new std::vector<float>);
0291     layerClustersMask->resize(clusterHandle->size(), 1.0);
0292     evt.put(std::move(layerClustersMask), "InitialLayerClustersMask");
0293   }
0294
0295   auto timeCl = std::make_unique<edm::ValueMap<std::pair<float, float>>>();
0296   edm::ValueMap<std::pair<float, float>>::Filler filler(*timeCl);
0297   filler.insert(clusterHandle, times.begin(), times.end());
0298   filler.fill();
0299   evt.put(std::move(timeCl), timeClname_);
0300
0301   algo_->reset();
0302 }
0303
0304 void HGCalLayerClusterProducer::setAlgoId() {
0305   if (detector_ == "EE") {
0306     algoId_ = reco::CaloCluster::hgcal_em;
0307   } else if (detector_ == "FH") {
0308     algoId_ = reco::CaloCluster::hgcal_had;
0309   } else if (detector_ == "BH") {
0310     algoId_ = reco::CaloCluster::hgcal_scintillator;
0311   } else if (detector_ == "HFNose") {
0312     algoId_ = reco::CaloCluster::hfnose;
0313   }
0314 }
0315
0316 #include "FWCore/Framework/interface/MakerMacros.h"
0317 DEFINE_FWK_MODULE(HGCalLayerClusterProducer);