Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-06-22 02:24:04

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