Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2025-01-14 23:17:10

0001 #include <Eigen/Core>
0002 #include <Eigen/Dense>
0003 #include <algorithm>
0004 #include <array>
0005 #include <cmath>
0006 #include <iostream>
0007 #include <memory>
0008 #include <string>
0009 #include <type_traits>
0010 #include <unordered_map>
0011 #include <utility>
0012 #include <vector>
0013 
0014 #include "FWCore/Framework/interface/stream/EDProducer.h"
0015 #include "FWCore/Utilities/interface/StreamID.h"
0016 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0017 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0018 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0019 #include "FWCore/Utilities/interface/EDGetToken.h"
0020 #include "FWCore/Utilities/interface/InputTag.h"
0021 #include "FWCore/Framework/interface/Event.h"
0022 #include "FWCore/Framework/interface/EventSetup.h"
0023 #include "FWCore/Utilities/interface/EDPutToken.h"
0024 
0025 #include "CondFormats/DataRecord/interface/HcalPFCutsRcd.h"
0026 #include "CondTools/Hcal/interface/HcalPFCutsHandler.h"
0027 
0028 #include "DataFormats/Common/interface/Handle.h"
0029 #include "DataFormats/DetId/interface/DetId.h"
0030 #include "DataFormats/ParticleFlowReco/interface/PFRecHitHostCollection.h"
0031 #include "DataFormats/ParticleFlowReco/interface/PFClusterHostCollection.h"
0032 #include "DataFormats/ParticleFlowReco/interface/PFRecHitFractionHostCollection.h"
0033 #include "RecoParticleFlow/PFClusterProducer/interface/PFCPositionCalculatorBase.h"
0034 
0035 class LegacyPFClusterProducer : public edm::stream::EDProducer<> {
0036 public:
0037   LegacyPFClusterProducer(edm::ParameterSet const& config)
0038       : pfClusterSoAToken_(consumes(config.getParameter<edm::InputTag>("src"))),
0039         pfRecHitFractionSoAToken_(consumes(config.getParameter<edm::InputTag>("src"))),
0040         InputPFRecHitSoA_Token_{consumes(config.getParameter<edm::InputTag>("PFRecHitsLabelIn"))},
0041         legacyPfClustersToken_(produces()),
0042         recHitsLabel_(consumes(config.getParameter<edm::InputTag>("recHitsSource"))),
0043         hcalCutsToken_(esConsumes<HcalPFCuts, HcalPFCutsRcd>(edm::ESInputTag("", "withTopo"))),
0044         cutsFromDB_(config.getParameter<bool>("usePFThresholdsFromDB")) {
0045     edm::ConsumesCollector cc = consumesCollector();
0046 
0047     //setup pf cluster builder if requested
0048     const edm::ParameterSet& pfcConf = config.getParameterSet("pfClusterBuilder");
0049     if (!pfcConf.empty()) {
0050       const auto& acConf = pfcConf.getParameterSet("positionCalc");
0051       if (!acConf.empty()) {
0052         const std::string& algoac = acConf.getParameter<std::string>("algoName");
0053         if (!algoac.empty())
0054           positionCalc_ = PFCPositionCalculatorFactory::get()->create(algoac, acConf, cc);
0055       }
0056 
0057       const auto& acConf2 = pfcConf.getParameterSet("allCellsPositionCalc");
0058       if (!acConf2.empty()) {
0059         const std::string& algoac = acConf2.getParameter<std::string>("algoName");
0060         if (!algoac.empty())
0061           allCellsPositionCalc_ = PFCPositionCalculatorFactory::get()->create(algoac, acConf2, cc);
0062       }
0063     }
0064   }
0065 
0066   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0067     edm::ParameterSetDescription desc;
0068     desc.add<edm::InputTag>("src", edm::InputTag("pfClusterSoAProducer"));
0069     desc.add<edm::InputTag>("PFRecHitsLabelIn", edm::InputTag("pfRecHitSoAProducerHCAL"));
0070     desc.add<edm::InputTag>("recHitsSource", edm::InputTag("legacyPFRecHitProducer"));
0071     desc.add<bool>("usePFThresholdsFromDB", true);
0072     {
0073       edm::ParameterSetDescription pfClusterBuilder;
0074       pfClusterBuilder.add<unsigned int>("maxIterations", 5);
0075       pfClusterBuilder.add<double>("minFracTot", 1e-20);
0076       pfClusterBuilder.add<double>("minFractionToKeep", 1e-7);
0077       pfClusterBuilder.add<bool>("excludeOtherSeeds", true);
0078       pfClusterBuilder.add<double>("showerSigma", 10.);
0079       pfClusterBuilder.add<double>("stoppingTolerance", 1e-8);
0080       pfClusterBuilder.add<double>("timeSigmaEB", 10.);
0081       pfClusterBuilder.add<double>("timeSigmaEE", 10.);
0082       pfClusterBuilder.add<double>("maxNSigmaTime", 10.);
0083       pfClusterBuilder.add<double>("minChi2Prob", 0.);
0084       pfClusterBuilder.add<bool>("clusterTimeResFromSeed", false);
0085       pfClusterBuilder.add<std::string>("algoName", "");
0086       pfClusterBuilder.add<edm::ParameterSetDescription>("positionCalcForConvergence", {});
0087       {
0088         edm::ParameterSetDescription validator;
0089         validator.add<std::string>("detector", "");
0090         validator.add<std::vector<int>>("depths", {});
0091         validator.add<std::vector<double>>("recHitEnergyNorm", {});
0092         std::vector<edm::ParameterSet> vDefaults(2);
0093         vDefaults[0].addParameter<std::string>("detector", "HCAL_BARREL1");
0094         vDefaults[0].addParameter<std::vector<int>>("depths", {1, 2, 3, 4});
0095         vDefaults[0].addParameter<std::vector<double>>("recHitEnergyNorm", {0.1, 0.2, 0.3, 0.3});
0096         vDefaults[1].addParameter<std::string>("detector", "HCAL_ENDCAP");
0097         vDefaults[1].addParameter<std::vector<int>>("depths", {1, 2, 3, 4, 5, 6, 7});
0098         vDefaults[1].addParameter<std::vector<double>>("recHitEnergyNorm", {0.1, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2});
0099         pfClusterBuilder.addVPSet("recHitEnergyNorms", validator, vDefaults);
0100       }
0101       {
0102         edm::ParameterSetDescription bar;
0103         bar.add<std::string>("algoName", "Basic2DGenericPFlowPositionCalc");
0104         bar.add<double>("minFractionInCalc", 1e-9);
0105         bar.add<int>("posCalcNCrystals", 5);
0106         {
0107           edm::ParameterSetDescription validator;
0108           validator.add<std::string>("detector", "");
0109           validator.add<std::vector<int>>("depths", {});
0110           validator.add<std::vector<double>>("logWeightDenominator", {});
0111           std::vector<edm::ParameterSet> vDefaults(2);
0112           vDefaults[0].addParameter<std::string>("detector", "HCAL_BARREL1");
0113           vDefaults[0].addParameter<std::vector<int>>("depths", {1, 2, 3, 4});
0114           vDefaults[0].addParameter<std::vector<double>>("logWeightDenominator", {0.1, 0.2, 0.3, 0.3});
0115           vDefaults[1].addParameter<std::string>("detector", "HCAL_ENDCAP");
0116           vDefaults[1].addParameter<std::vector<int>>("depths", {1, 2, 3, 4, 5, 6, 7});
0117           vDefaults[1].addParameter<std::vector<double>>("logWeightDenominator", {0.1, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2});
0118           bar.addVPSet("logWeightDenominatorByDetector", validator, vDefaults);
0119         }
0120         bar.add<double>("minAllowedNormalization", 1e-9);
0121         bar.add<edm::ParameterSetDescription>("timeResolutionCalcBarrel", {});
0122         bar.add<edm::ParameterSetDescription>("timeResolutionCalcEndcap", {});
0123         pfClusterBuilder.add("positionCalc", bar);
0124       }
0125       {
0126         edm::ParameterSetDescription bar;
0127         bar.add<std::string>("algoName", "Basic2DGenericPFlowPositionCalc");
0128         bar.add<double>("minFractionInCalc", 1e-9);
0129         bar.add<int>("posCalcNCrystals", -1);
0130         {
0131           edm::ParameterSetDescription validator;
0132           validator.add<std::string>("detector", "");
0133           validator.add<std::vector<int>>("depths", {});
0134           validator.add<std::vector<double>>("logWeightDenominator", {});
0135           std::vector<edm::ParameterSet> vDefaults(2);
0136           vDefaults[0].addParameter<std::string>("detector", "HCAL_BARREL1");
0137           vDefaults[0].addParameter<std::vector<int>>("depths", {1, 2, 3, 4});
0138           vDefaults[0].addParameter<std::vector<double>>("logWeightDenominator", {0.1, 0.2, 0.3, 0.3});
0139           vDefaults[1].addParameter<std::string>("detector", "HCAL_ENDCAP");
0140           vDefaults[1].addParameter<std::vector<int>>("depths", {1, 2, 3, 4, 5, 6, 7});
0141           vDefaults[1].addParameter<std::vector<double>>("logWeightDenominator", {0.1, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2});
0142           bar.addVPSet("logWeightDenominatorByDetector", validator, vDefaults);
0143         }
0144         bar.add<double>("minAllowedNormalization", 1e-9);
0145         bar.add<edm::ParameterSetDescription>("timeResolutionCalcBarrel", {});
0146         bar.add<edm::ParameterSetDescription>("timeResolutionCalcEndcap", {});
0147         pfClusterBuilder.add("allCellsPositionCalc", bar);
0148       }
0149       {
0150         edm::ParameterSetDescription bar;
0151         bar.add<double>("corrTermLowE", 0.);
0152         bar.add<double>("threshLowE", 6.);
0153         bar.add<double>("noiseTerm", 21.86);
0154         bar.add<double>("constantTermLowE", 4.24);
0155         bar.add<double>("noiseTermLowE", 8.);
0156         bar.add<double>("threshHighE", 15.);
0157         bar.add<double>("constantTerm", 2.82);
0158         pfClusterBuilder.add("timeResolutionCalcBarrel", bar);
0159       }
0160       {
0161         edm::ParameterSetDescription bar;
0162         bar.add<double>("corrTermLowE", 0.);
0163         bar.add<double>("threshLowE", 6.);
0164         bar.add<double>("noiseTerm", 21.86);
0165         bar.add<double>("constantTermLowE", 4.24);
0166         bar.add<double>("noiseTermLowE", 8.);
0167         bar.add<double>("threshHighE", 15.);
0168         bar.add<double>("constantTerm", 2.82);
0169         pfClusterBuilder.add("timeResolutionCalcEndcap", bar);
0170       }
0171       {
0172         edm::ParameterSetDescription bar;
0173         pfClusterBuilder.add("positionReCalc", bar);
0174       }
0175       {
0176         edm::ParameterSetDescription bar;
0177         pfClusterBuilder.add("energyCorrector", bar);
0178       }
0179       desc.add("pfClusterBuilder", pfClusterBuilder);
0180     }
0181     descriptions.addWithDefaultLabel(desc);
0182   }
0183 
0184 private:
0185   void produce(edm::Event&, const edm::EventSetup&) override;
0186   const edm::EDGetTokenT<reco::PFClusterHostCollection> pfClusterSoAToken_;
0187   const edm::EDGetTokenT<reco::PFRecHitFractionHostCollection> pfRecHitFractionSoAToken_;
0188   const edm::EDGetTokenT<reco::PFRecHitHostCollection> InputPFRecHitSoA_Token_;
0189   const edm::EDPutTokenT<reco::PFClusterCollection> legacyPfClustersToken_;
0190   const edm::EDGetTokenT<reco::PFRecHitCollection> recHitsLabel_;
0191   const edm::ESGetToken<HcalPFCuts, HcalPFCutsRcd> hcalCutsToken_;
0192   const bool cutsFromDB_;
0193   // the actual algorithm
0194   std::unique_ptr<PFCPositionCalculatorBase> positionCalc_;
0195   std::unique_ptr<PFCPositionCalculatorBase> allCellsPositionCalc_;
0196 };
0197 
0198 void LegacyPFClusterProducer::produce(edm::Event& event, const edm::EventSetup& setup) {
0199   const reco::PFRecHitHostCollection& pfRecHits = event.get(InputPFRecHitSoA_Token_);
0200 
0201   HcalPFCuts const* paramPF = cutsFromDB_ ? &setup.getData(hcalCutsToken_) : nullptr;
0202 
0203   auto const& pfClusterSoA = event.get(pfClusterSoAToken_).const_view();
0204   auto const& pfRecHitFractionSoA = event.get(pfRecHitFractionSoAToken_).const_view();
0205 
0206   int nRH = 0;
0207   if (pfRecHits->metadata().size() != 0)
0208     nRH = pfRecHits.view().size();
0209   reco::PFClusterCollection out;
0210   out.reserve(nRH);
0211 
0212   auto const rechitsHandle = event.getHandle(recHitsLabel_);
0213 
0214   if (nRH != 0) {
0215     // Build PFClusters in legacy format
0216     std::vector<int> nTopoSeeds(nRH, 0);
0217 
0218     for (int i = 0; i < pfClusterSoA.nSeeds(); i++) {
0219       nTopoSeeds[pfClusterSoA[i].topoId()]++;
0220     }
0221 
0222     // Looping over SoA PFClusters to produce legacy PFCluster collection
0223     for (int i = 0; i < pfClusterSoA.nSeeds(); i++) {
0224       unsigned int n = pfClusterSoA[i].seedRHIdx();
0225       reco::PFCluster temp;
0226       temp.setSeed((*rechitsHandle)[n].detId());  // Pulling the detId of this PFRecHit from the legacy format input
0227       int offset = pfClusterSoA[i].rhfracOffset();
0228       for (int k = offset; k < (offset + pfClusterSoA[i].rhfracSize()) && k >= 0;
0229            k++) {  // Looping over PFRecHits in the same topo cluster
0230         if (pfRecHitFractionSoA[k].pfrhIdx() < nRH && pfRecHitFractionSoA[k].pfrhIdx() > -1 &&
0231             pfRecHitFractionSoA[k].frac() > 0.0) {
0232           const reco::PFRecHitRef& refhit = reco::PFRecHitRef(rechitsHandle, pfRecHitFractionSoA[k].pfrhIdx());
0233           temp.addRecHitFraction(reco::PFRecHitFraction(refhit, pfRecHitFractionSoA[k].frac()));
0234         }
0235       }
0236 
0237       // Now PFRecHitFraction of this PFCluster is set. Now compute calculateAndSetPosition (energy, position etc)
0238       if (nTopoSeeds[pfClusterSoA[i].topoId()] == 1 && allCellsPositionCalc_) {
0239         allCellsPositionCalc_->calculateAndSetPosition(temp, paramPF);
0240       } else {
0241         positionCalc_->calculateAndSetPosition(temp, paramPF);
0242       }
0243       out.emplace_back(std::move(temp));
0244     }
0245   }
0246 
0247   event.emplace(legacyPfClustersToken_, std::move(out));
0248 }
0249 
0250 #include "FWCore/Framework/interface/MakerMacros.h"
0251 DEFINE_FWK_MODULE(LegacyPFClusterProducer);