File indexing completed on 2023-10-25 10:00:11
0001
0002
0003
0004
0005
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
0039
0040
0041
0042
0043
0044 HGCalLayerClusterProducer(const edm::ParameterSet&);
0045 ~HGCalLayerClusterProducer() override {}
0046
0047
0048
0049
0050
0051 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0052
0053
0054
0055
0056
0057
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
0073 std::vector<double> thresholdW0_;
0074 double positionDeltaRho2_;
0075 hgcal::RecHitTools rhtools_;
0076 edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeomToken_;
0077
0078
0079
0080
0081 void setAlgoId();
0082
0083
0084
0085
0086
0087
0088
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
0095
0096
0097
0098
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")),
0108 timeClname_(ps.getParameter<std::string>("timeClname")),
0109 hitsTime_(ps.getParameter<unsigned int>("nHitsTime")),
0110 caloGeomToken_(consumesCollector().esConsumes<CaloGeometry, CaloGeometryRecord>()) {
0111 setAlgoId();
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
0128 produces<edm::ValueMap<std::pair<float, float>>>(timeClname_);
0129 }
0130
0131 void HGCalLayerClusterProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0132
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
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
0168 const HGCRecHit* rechit = hitmap[hit.first];
0169
0170 const GlobalPoint position(rhtools_.getPosition(rechit->detid()));
0171
0172 if (thick != -1) {
0173
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 {
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
0211 const HGCRecHit* rechit = hitmap[hit.first];
0212
0213 float rhTimeE = rechit->timeError();
0214
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
0235
0236
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 {
0284 algoId_ = reco::CaloCluster::hgcal_had;
0285 }
0286 }
0287
0288 #include "FWCore/Framework/interface/MakerMacros.h"
0289 DEFINE_FWK_MODULE(HGCalLayerClusterProducer);