HGCalRecProducers/plugins/HGCalLayerClusterProducer.cc

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