PFClusterProducer/plugins/PFClusterProducer.cc

0001 #include "FWCore/Framework/interface/Event.h"
0002 #include "FWCore/Framework/interface/stream/EDProducer.h"
0003 #include "RecoParticleFlow/PFClusterProducer/interface/InitialClusteringStepBase.h"
0004 #include "RecoParticleFlow/PFClusterProducer/interface/PFCPositionCalculatorBase.h"
0005 #include "RecoParticleFlow/PFClusterProducer/interface/PFClusterBuilderBase.h"
0006 #include "RecoParticleFlow/PFClusterProducer/interface/PFClusterEnergyCorrectorBase.h"
0007 #include "RecoParticleFlow/PFClusterProducer/interface/RecHitTopologicalCleanerBase.h"
0008 #include "RecoParticleFlow/PFClusterProducer/interface/SeedFinderBase.h"
0009 #include "CondFormats/DataRecord/interface/HcalPFCutsRcd.h"
0010 #include "CondTools/Hcal/interface/HcalPFCutsHandler.h"
0011 #include <memory>
0012
0013 class PFClusterProducer : public edm::stream::EDProducer<> {
0014   typedef RecHitTopologicalCleanerBase RHCB;
0015   typedef InitialClusteringStepBase ICSB;
0016   typedef PFClusterBuilderBase PFCBB;
0017   typedef PFCPositionCalculatorBase PosCalc;
0018
0019 public:
0020   PFClusterProducer(const edm::ParameterSet&);
0021   ~PFClusterProducer() override = default;
0022
0023   void beginRun(const edm::Run&, const edm::EventSetup&) override;
0024   void produce(edm::Event&, const edm::EventSetup&) override;
0025   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0026
0027 private:
0028   // inputs
0029   edm::EDGetTokenT<reco::PFRecHitCollection> _rechitsLabel;
0030   edm::ESGetToken<HcalPFCuts, HcalPFCutsRcd> hcalCutsToken_;
0031   bool cutsFromDB;
0032   HcalPFCuts const* paramPF = nullptr;
0033
0034   // options
0035   const bool _prodInitClusters;
0036   // the actual algorithm
0037   std::vector<std::unique_ptr<RecHitTopologicalCleanerBase>> _cleaners;
0038   std::vector<std::unique_ptr<RecHitTopologicalCleanerBase>> _seedcleaners;
0039   std::unique_ptr<SeedFinderBase> _seedFinder;
0040   std::unique_ptr<InitialClusteringStepBase> _initialClustering;
0041   std::unique_ptr<PFClusterBuilderBase> _pfClusterBuilder;
0042   std::unique_ptr<PFCPositionCalculatorBase> _positionReCalc;
0043   std::unique_ptr<PFClusterEnergyCorrectorBase> _energyCorrector;
0044 };
0045
0046 #include "FWCore/Framework/interface/MakerMacros.h"
0047 DEFINE_FWK_MODULE(PFClusterProducer);
0048
0049 void PFClusterProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0050   edm::ParameterSetDescription desc;
0051   desc.add<edm::InputTag>("recHitsSource", {});
0052   desc.add<bool>("usePFThresholdsFromDB", false);
0053   {
0054     edm::ParameterSetDescription psd;
0055     psd.add<std::string>("algoName", "");
0056     desc.addVPSet("recHitCleaners", psd, {});
0057   }
0058   {
0059     edm::ParameterSetDescription psd;
0060     psd.add<std::string>("algoName", "");
0061     psd.add<std::vector<std::string>>("RecHitFlagsToBeExcluded", {});
0062     desc.addVPSet("seedCleaners", psd, {});
0063   }
0064   {
0065     edm::ParameterSetDescription pset;
0066     pset.add<std::string>("algoName", "");
0067     pset.add<int>("nNeighbours", 0);
0068     {
0069       edm::ParameterSetDescription psd;
0070       psd.add<std::string>("detector", "");
0071       psd.addNode((edm::ParameterDescription<double>("seedingThreshold", 0, true) and
0072                    edm::ParameterDescription<double>("seedingThresholdPt", 0, true)) xor
0073                   (edm::ParameterDescription<std::vector<int>>("depths", {}, true) and
0074                    edm::ParameterDescription<std::vector<double>>("seedingThreshold", {}, true) and
0075                    edm::ParameterDescription<std::vector<double>>("seedingThresholdPt", {}, true)));
0076       pset.addVPSet("thresholdsByDetector", psd, {});
0077     }
0078     desc.add<edm::ParameterSetDescription>("seedFinder", pset);
0079   }
0080   {
0081     edm::ParameterSetDescription pset;
0082     pset.add<std::string>("algoName", "");
0083     {
0084       edm::ParameterSetDescription psd;
0085       psd.add<std::string>("detector", "");
0086       psd.addNode((edm::ParameterDescription<double>("gatheringThreshold", 0, true) and
0087                    edm::ParameterDescription<double>("gatheringThresholdPt", 0, true)) xor
0088                   (edm::ParameterDescription<std::vector<int>>("depths", {}, true) and
0089                    edm::ParameterDescription<std::vector<double>>("gatheringThreshold", {}, true) and
0090                    edm::ParameterDescription<std::vector<double>>("gatheringThresholdPt", {}, true)));
0091       pset.addVPSet("thresholdsByDetector", psd, {});
0092     }
0093     pset.add<bool>("useCornerCells", false);
0094     pset.add<edm::InputTag>("clusterSrc", {});
0095     pset.add<bool>("filterByTracksterIteration", false);
0096     pset.add<bool>("filterByTracksterPID", false);
0097     pset.add<std::vector<int>>("filter_on_categories", {});
0098     pset.add<std::vector<int>>("filter_on_iterations", {});
0099     pset.add<double>("pid_threshold", 0);
0100     pset.add<edm::InputTag>("tracksterSrc", {});
0101     pset.add<double>("exclusiveFraction", 0);
0102     pset.add<double>("invisibleFraction", 0);
0103     pset.add<bool>("maxDistanceFilter", false);
0104     pset.add<double>("maxDistance", 0);
0105     pset.add<double>("maxDforTimingSquared", 0);
0106     pset.add<double>("timeOffset", 0);
0107     pset.add<uint32_t>("minNHitsforTiming", 0);
0108     pset.add<bool>("useMCFractionsForExclEnergy", false);
0109     pset.add<std::vector<double>>("hadronCalib", {});
0110     pset.add<std::vector<double>>("egammaCalib", {});
0111     pset.add<double>("calibMinEta", 0);
0112     pset.add<double>("calibMaxEta", 0);
0113     pset.add<edm::InputTag>("simClusterSrc", {});
0114     desc.add<edm::ParameterSetDescription>("initialClusteringStep", pset);
0115   }
0116   {
0117     edm::ParameterSetDescription pset;
0118     pset.add<std::string>("algoName", "");
0119     {
0120       edm::ParameterSetDescription pset2;
0121       pset2.add<std::string>("algoName", "");
0122       pset2.add<double>("minFractionInCalc", 0);
0123       pset2.add<int>("posCalcNCrystals", -1);
0124       {
0125         edm::ParameterSetDescription psd;
0126         psd.add<std::string>("detector", "");
0127         psd.add<std::vector<int>>("depths", {});
0128         psd.add<std::vector<double>>("logWeightDenominator", {});
0129         pset2.addVPSet("logWeightDenominatorByDetector", psd, {});
0130       }
0131       pset2.add<double>("logWeightDenominator", 0);
0132       pset2.add<double>("minAllowedNormalization", 0);
0133       {
0134         edm::ParameterSetDescription pset3;
0135         pset3.add<double>("constantTerm", 0);
0136         pset3.add<double>("constantTermLowE", 0);
0137         pset3.add<double>("corrTermLowE", 0);
0138         pset3.add<double>("noiseTerm", 0);
0139         pset3.add<double>("noiseTermLowE", 0);
0140         pset3.add<double>("threshHighE", -1.);
0141         pset3.add<double>("threshLowE", -1.);
0142         pset2.add<edm::ParameterSetDescription>("timeResolutionCalcBarrel", pset3);
0143         pset2.add<edm::ParameterSetDescription>("timeResolutionCalcEndcap", pset3);
0144       }
0145       pset.add<edm::ParameterSetDescription>("allCellsPositionCalc", pset2);
0146       pset.add<edm::ParameterSetDescription>("positionCalc", pset2);
0147     }
0148     pset.add<double>("minFractionToKeep", 0);
0149     pset.add<double>("nSigmaEta", 0);
0150     pset.add<double>("nSigmaPhi", 0);
0151     pset.add<bool>("excludeOtherSeeds", false);
0152     pset.add<uint32_t>("maxIterations", 0);
0153     pset.add<double>("minFracTot", 0);
0154     {
0155       edm::ParameterSetDescription pset2;
0156       pset2.add<std::string>("algoName", "");
0157       pset2.add<double>("minFractionInCalc", 0);
0158       pset2.add<double>("T0_EB", 0);
0159       pset2.add<double>("T0_EE", 0);
0160       pset2.add<double>("T0_ES", 0);
0161       pset2.add<double>("W0", 0);
0162       pset2.add<double>("X0", 0);
0163       pset2.add<double>("minAllowedNormalization", 0);
0164       pset2.add<edm::ParameterSetDescription>("timeResolutionCalc", {});
0165       pset.add<edm::ParameterSetDescription>("positionCalcForConvergence", pset2);
0166     }
0167     {
0168       edm::ParameterSetDescription psd;
0169       psd.add<std::string>("detector", "");
0170       psd.addNode(edm::ParameterDescription<double>("recHitEnergyNorm", 0, true) xor
0171                   (edm::ParameterDescription<std::vector<int>>("depths", {}, true) and
0172                    edm::ParameterDescription<std::vector<double>>("recHitEnergyNorm", {}, true)));
0173       pset.addVPSet("recHitEnergyNorms", psd, {});
0174     }
0175     pset.add<double>("showerSigma", 1.5);
0176     pset.add<double>("stoppingTolerance", 1e-08);
0177     pset.add<bool>("clusterTimeResFromSeed", false);
0178     pset.add<double>("maxNSigmaTime", 10.0);
0179     pset.add<double>("minChi2Prob", 0);
0180     {
0181       edm::ParameterSetDescription pset2;
0182       pset2.add<double>("constantTerm", 0);
0183       pset2.add<double>("constantTermLowE", 0);
0184       pset2.add<double>("corrTermLowE", 0);
0185       pset2.add<double>("noiseTerm", 0);
0186       pset2.add<double>("noiseTermLowE", 0);
0187       pset2.add<double>("threshHighE", -1.);
0188       pset2.add<double>("threshLowE", -1.);
0189       pset.add<edm::ParameterSetDescription>("timeResolutionCalcBarrel", pset2);
0190       pset.add<edm::ParameterSetDescription>("timeResolutionCalcEndcap", pset2);
0191     }
0192     pset.add<double>("timeSigmaEB", 10.0);
0193     pset.add<double>("timeSigmaEE", 10.0);
0194     desc.add<edm::ParameterSetDescription>("pfClusterBuilder", pset);
0195   }
0196   {
0197     edm::ParameterSetDescription pset;
0198     pset.add<std::string>("algoName", "");
0199     pset.add<double>("minFractionInCalc", 0);
0200     pset.add<bool>("updateTiming", false);
0201     pset.add<double>("T0_EB", 0);
0202     pset.add<double>("T0_EE", 0);
0203     pset.add<double>("T0_ES", 0);
0204     pset.add<double>("W0", 0);
0205     pset.add<double>("X0", 0);
0206     pset.add<double>("minAllowedNormalization", 0);
0207     pset.add<edm::ParameterSetDescription>("timeResolutionCalc", {});
0208     desc.add<edm::ParameterSetDescription>("positionReCalc", pset);
0209   }
0210   desc.add<edm::ParameterSetDescription>("energyCorrector", {});
0211   descriptions.addWithDefaultLabel(desc);
0212 }
0213
0214 #ifdef PFLOW_DEBUG
0215 #define LOGVERB(x) edm::LogVerbatim(x)
0216 #define LOGWARN(x) edm::LogWarning(x)
0217 #define LOGERR(x) edm::LogError(x)
0218 #define LOGDRESSED(x) edm::LogInfo(x)
0219 #else
0220 #define LOGVERB(x) LogTrace(x)
0221 #define LOGWARN(x) edm::LogWarning(x)
0222 #define LOGERR(x) edm::LogError(x)
0223 #define LOGDRESSED(x) LogDebug(x)
0224 #endif
0225
0226 PFClusterProducer::PFClusterProducer(const edm::ParameterSet& conf)
0227     : _prodInitClusters(conf.getUntrackedParameter<bool>("prodInitialClusters", false)) {
0228   _rechitsLabel = consumes<reco::PFRecHitCollection>(conf.getParameter<edm::InputTag>("recHitsSource"));
0229   cutsFromDB = conf.getParameter<bool>("usePFThresholdsFromDB");
0230   edm::ConsumesCollector cc = consumesCollector();
0231
0232   if (cutsFromDB) {
0233     hcalCutsToken_ = esConsumes<HcalPFCuts, HcalPFCutsRcd, edm::Transition::BeginRun>(edm::ESInputTag("", "withTopo"));
0234   }
0235
0236   //setup rechit cleaners
0237   const edm::VParameterSet& cleanerConfs = conf.getParameterSetVector("recHitCleaners");
0238   for (const auto& conf : cleanerConfs) {
0239     const std::string& cleanerName = conf.getParameter<std::string>("algoName");
0240     _cleaners.emplace_back(RecHitTopologicalCleanerFactory::get()->create(cleanerName, conf, cc));
0241   }
0242
0243   const auto& seedcleanerConfs = conf.getParameterSetVector("seedCleaners");
0244   if (!seedcleanerConfs.empty()) {
0245     for (const auto& conf : seedcleanerConfs) {
0246       const std::string& seedcleanerName = conf.getParameter<std::string>("algoName");
0247       _seedcleaners.emplace_back(RecHitTopologicalCleanerFactory::get()->create(seedcleanerName, conf, cc));
0248     }
0249   }
0250
0251   // setup seed finding
0252   const edm::ParameterSet& sfConf = conf.getParameterSet("seedFinder");
0253   const std::string& sfName = sfConf.getParameter<std::string>("algoName");
0254   _seedFinder = SeedFinderFactory::get()->create(sfName, sfConf);
0255   //setup topo cluster builder
0256   const edm::ParameterSet& initConf = conf.getParameterSet("initialClusteringStep");
0257   const std::string& initName = initConf.getParameter<std::string>("algoName");
0258   _initialClustering = InitialClusteringStepFactory::get()->create(initName, initConf, cc);
0259   //setup pf cluster builder if requested
0260   const edm::ParameterSet& pfcConf = conf.getParameterSet("pfClusterBuilder");
0261   if (!pfcConf.empty()) {
0262     const std::string& pfcName = pfcConf.getParameter<std::string>("algoName");
0263     if (!pfcName.empty())
0264       _pfClusterBuilder = PFClusterBuilderFactory::get()->create(pfcName, pfcConf, cc);
0265   }
0266   //setup (possible) recalcuation of positions
0267   const edm::ParameterSet& pConf = conf.getParameterSet("positionReCalc");
0268   if (!pConf.empty()) {
0269     const std::string& pName = pConf.getParameter<std::string>("algoName");
0270     if (!pName.empty())
0271       _positionReCalc = PFCPositionCalculatorFactory::get()->create(pName, pConf, cc);
0272   }
0273   // see if new need to apply corrections, setup if there.
0274   const edm::ParameterSet& cConf = conf.getParameterSet("energyCorrector");
0275   if (!cConf.empty()) {
0276     const std::string& cName = cConf.getParameter<std::string>("algoName");
0277     if (!cName.empty())
0278       _energyCorrector = PFClusterEnergyCorrectorFactory::get()->create(cName, cConf);
0279   }
0280
0281   if (_prodInitClusters) {
0282     produces<reco::PFClusterCollection>("initialClusters");
0283   }
0284   produces<reco::PFClusterCollection>();
0285 }
0286
0287 void PFClusterProducer::beginRun(const edm::Run& run, const edm::EventSetup& es) {
0288   if (cutsFromDB) {
0289     paramPF = &es.getData(hcalCutsToken_);
0290   }
0291   _initialClustering->update(es);
0292   if (_pfClusterBuilder)
0293     _pfClusterBuilder->update(es);
0294   if (_positionReCalc)
0295     _positionReCalc->update(es);
0296   for (const auto& cleaner : _cleaners)
0297     cleaner->update(es);
0298   for (const auto& cleaner : _seedcleaners)
0299     cleaner->update(es);
0300 }
0301
0302 void PFClusterProducer::produce(edm::Event& e, const edm::EventSetup& es) {
0303   _initialClustering->reset();
0304   if (_pfClusterBuilder)
0305     _pfClusterBuilder->reset();
0306
0307   edm::Handle<reco::PFRecHitCollection> rechits;
0308   e.getByToken(_rechitsLabel, rechits);
0309
0310   _initialClustering->updateEvent(e);
0311
0312   std::vector<bool> mask(rechits->size(), true);
0313   for (const auto& cleaner : _cleaners) {
0314     cleaner->clean(rechits, mask);
0315   }
0316
0317   // no seeding on these hits
0318   std::vector<bool> seedmask = mask;
0319   for (const auto& cleaner : _seedcleaners) {
0320     cleaner->clean(rechits, seedmask);
0321   }
0322
0323   std::vector<bool> seedable(rechits->size(), false);
0324   _seedFinder->findSeeds(rechits, seedmask, seedable, paramPF);
0325
0326   auto initialClusters = std::make_unique<reco::PFClusterCollection>();
0327   _initialClustering->buildClusters(rechits, mask, seedable, *initialClusters, paramPF);
0328   LOGVERB("PFClusterProducer::produce()") << *_initialClustering;
0329
0330   auto pfClusters = std::make_unique<reco::PFClusterCollection>();
0331   pfClusters = std::make_unique<reco::PFClusterCollection>();
0332   if (_pfClusterBuilder) {  // if we've defined a re-clustering step execute it
0333     _pfClusterBuilder->buildClusters(*initialClusters, seedable, *pfClusters, paramPF);
0334     LOGVERB("PFClusterProducer::produce()") << *_pfClusterBuilder;
0335   } else {
0336     pfClusters->insert(pfClusters->end(), initialClusters->begin(), initialClusters->end());
0337   }
0338
0339   if (_positionReCalc) {
0340     _positionReCalc->calculateAndSetPositions(*pfClusters, paramPF);
0341   }
0342
0343   if (_energyCorrector) {
0344     _energyCorrector->correctEnergies(*pfClusters);
0345   }
0346
0347   if (_prodInitClusters)
0348     e.put(std::move(initialClusters), "initialClusters");
0349   e.put(std::move(pfClusters));
0350 }