PFClusterProducer/plugins/CorrectedECALPFClusterProducer.cc

0001 #ifndef __CorrectedECALPFClusterProducer__
0002 #define __CorrectedECALPFClusterProducer__
0003
0004 // user include files
0005 #include <memory>
0006
0007 #include "FWCore/Framework/interface/Frameworkfwd.h"
0008 #include "FWCore/Framework/interface/stream/EDProducer.h"
0009
0010 #include "FWCore/Framework/interface/Event.h"
0011 #include "FWCore/Framework/interface/MakerMacros.h"
0012 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0013 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0014 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0015
0016 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
0017 #include "DataFormats/ParticleFlowReco/interface/PFClusterFwd.h"
0018
0019 #include "RecoParticleFlow/PFClusterProducer/interface/PFClusterEMEnergyCorrector.h"
0020 #include "RecoParticleFlow/PFClusterTools/interface/LinkByRecHit.h"
0021 #include "DataFormats/Math/interface/deltaPhi.h"
0022
0023 namespace {
0024   typedef reco::PFCluster::EEtoPSAssociation::value_type EEPSPair;
0025   bool sortByKey(const EEPSPair& a, const EEPSPair& b) { return a.first < b.first; }
0026
0027   // why, why -1 and not <double>::max() ????
0028   double testPreshowerDistance(reco::PFCluster const& eeclus, reco::PFCluster const& psclus) {
0029     auto const& pspos = psclus.positionREP();
0030     auto const& eepos = eeclus.positionREP();
0031     // lazy continue based on geometry
0032     if (eeclus.z() * psclus.z() < 0)
0033       return -1.0;
0034     auto deta = std::abs(eepos.eta() - pspos.eta());
0035     if (deta > 0.3)
0036       return -1.0;
0037     auto dphi = std::abs(deltaPhi(eepos.phi(), pspos.phi()));
0038     if (dphi > 0.6)
0039       return -1.0;
0040     return LinkByRecHit::testECALAndPSByRecHit(eeclus, psclus, false);
0041   }
0042 }  // namespace
0043
0044 class CorrectedECALPFClusterProducer : public edm::stream::EDProducer<> {
0045 public:
0046   CorrectedECALPFClusterProducer(const edm::ParameterSet& conf)
0047       : minimumPSEnergy_(conf.getParameter<double>("minimumPSEnergy")), skipPS_(conf.getParameter<bool>("skipPS")) {
0048     const edm::InputTag& inputECAL = conf.getParameter<edm::InputTag>("inputECAL");
0049     inputECAL_ = consumes<reco::PFClusterCollection>(inputECAL);
0050
0051     const edm::InputTag& inputPS = conf.getParameter<edm::InputTag>("inputPS");
0052     if (!skipPS_)
0053       inputPS_ = consumes<reco::PFClusterCollection>(inputPS);
0054
0055     const edm::ParameterSet& corConf = conf.getParameterSet("energyCorrector");
0056     corrector_ = std::make_unique<PFClusterEMEnergyCorrector>(corConf, consumesCollector());
0057
0058     produces<reco::PFCluster::EEtoPSAssociation>();
0059     produces<reco::PFClusterCollection>();
0060   }
0061
0062   void produce(edm::Event& e, const edm::EventSetup& es) override;
0063   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0064
0065 private:
0066   const double minimumPSEnergy_;
0067   const bool skipPS_;
0068   std::unique_ptr<PFClusterEMEnergyCorrector> corrector_;
0069   edm::EDGetTokenT<reco::PFClusterCollection> inputECAL_;
0070   edm::EDGetTokenT<reco::PFClusterCollection> inputPS_;
0071 };
0072
0073 DEFINE_FWK_MODULE(CorrectedECALPFClusterProducer);
0074
0075 void CorrectedECALPFClusterProducer::produce(edm::Event& e, const edm::EventSetup& es) {
0076   auto clusters_out = std::make_unique<reco::PFClusterCollection>();
0077   auto association_out = std::make_unique<reco::PFCluster::EEtoPSAssociation>();
0078
0079   edm::Handle<reco::PFClusterCollection> handleECAL;
0080   e.getByToken(inputECAL_, handleECAL);
0081   edm::Handle<reco::PFClusterCollection> handlePS;
0082   if (!skipPS_)
0083     e.getByToken(inputPS_, handlePS);
0084
0085   auto const& ecals = *handleECAL;
0086
0087   clusters_out->reserve(ecals.size());
0088   association_out->reserve(ecals.size());
0089   clusters_out->insert(clusters_out->end(), ecals.begin(), ecals.end());
0090
0091   //build the EE->PS association
0092   if (!skipPS_) {
0093     auto const& pss = *handlePS;
0094     for (unsigned i = 0; i < pss.size(); ++i) {
0095       switch (pss[i].layer()) {  // just in case this isn't the ES...
0096         case PFLayer::PS1:
0097         case PFLayer::PS2:
0098           break;
0099         default:
0100           continue;
0101       }
0102       if (pss[i].energy() < minimumPSEnergy_)
0103         continue;
0104       int eematch = -1;
0105       auto min_dist = std::numeric_limits<double>::max();
0106       for (size_t ic = 0; ic < ecals.size(); ++ic) {
0107         if (ecals[ic].layer() != PFLayer::ECAL_ENDCAP)
0108           continue;
0109         auto dist = testPreshowerDistance(ecals[ic], pss[i]);
0110         if (dist == -1.0)
0111           dist = std::numeric_limits<double>::max();
0112         if (dist < min_dist) {
0113           eematch = ic;
0114           min_dist = dist;
0115         }
0116       }  // loop on EE clusters
0117       if (eematch >= 0) {
0118         edm::Ptr<reco::PFCluster> psclus(handlePS, i);
0119         association_out->push_back(std::make_pair(eematch, psclus));
0120       }
0121     }
0122   }
0123   std::sort(association_out->begin(), association_out->end(), sortByKey);
0124
0125   corrector_->correctEnergies(e, es, *association_out, *clusters_out);
0126
0127   association_out->shrink_to_fit();
0128
0129   e.put(std::move(association_out));
0130   e.put(std::move(clusters_out));
0131 }
0132
0133 void CorrectedECALPFClusterProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0134   edm::ParameterSetDescription desc;
0135   desc.add<double>("minimumPSEnergy", 0.0);
0136   desc.ifValue(
0137       edm::ParameterDescription<bool>("skipPS", false, true),
0138       true >> (edm::ParameterDescription<edm::InputTag>("inputPS", edm::InputTag(""), true)) or
0139           false >> (edm::ParameterDescription<edm::InputTag>("inputPS", edm::InputTag("particleFlowClusterPS"), true)));
0140   {
0141     edm::ParameterSetDescription psd0;
0142     psd0.add<bool>("applyCrackCorrections", false);
0143     psd0.add<bool>("applyMVACorrections", false);
0144     psd0.add<bool>("srfAwareCorrection", false);
0145     psd0.add<bool>("setEnergyUncertainty", false);
0146     psd0.add<bool>("autoDetectBunchSpacing", true);
0147     psd0.add<int>("bunchSpacing", 25);
0148     psd0.add<double>("maxPtForMVAEvaluation", -99.);
0149     psd0.add<edm::InputTag>("recHitsEBLabel", edm::InputTag("ecalRecHit", "EcalRecHitsEB"));
0150     psd0.add<edm::InputTag>("recHitsEELabel", edm::InputTag("ecalRecHit", "EcalRecHitsEE"));
0151     psd0.add<edm::InputTag>("ebSrFlagLabel", edm::InputTag("ecalDigis"));
0152     psd0.add<edm::InputTag>("eeSrFlagLabel", edm::InputTag("ecalDigis"));
0153     desc.add<edm::ParameterSetDescription>("energyCorrector", psd0);
0154   }
0155   desc.add<edm::InputTag>("inputECAL", edm::InputTag("particleFlowClusterECALUncorrected"));
0156   descriptions.add("particleFlowClusterECAL", desc);
0157 }
0158
0159 #endif