plugins/SimMappers/RealisticSimClusterMapper.cc

0001 /////////////////////////
0002 // Author: Felice Pantaleo
0003 // Date:   30/06/2017
0004 // Email: felice@cern.ch
0005 /////////////////////////
0006
0007 #include "RealisticCluster.h"
0008
0009 #include "DataFormats/ForwardDetId/interface/HGCalDetId.h"
0010 #include "DataFormats/ParticleFlowReco/interface/PFRecHitFraction.h"
0011 #include "FWCore/Framework/interface/Event.h"
0012 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0013 #include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h"
0014 #include "RecoLocalCalo/HGCalRecProducers/interface/ComputeClusterTime.h"
0015 #include "RecoParticleFlow/PFClusterProducer/interface/InitialClusteringStepBase.h"
0016 #include "RecoParticleFlow/PFClusterProducer/plugins/SimMappers/RealisticHitToClusterAssociator.h"
0017 #include "SimDataFormats/CaloAnalysis/interface/SimCluster.h"
0018 #include "SimDataFormats/CaloAnalysis/interface/SimClusterFwd.h"
0019
0020 #include <unordered_map>
0021
0022 class RealisticSimClusterMapper : public InitialClusteringStepBase {
0023 public:
0024   RealisticSimClusterMapper(const edm::ParameterSet& conf, edm::ConsumesCollector& cc)
0025       : InitialClusteringStepBase(conf, cc),
0026         invisibleFraction_(conf.getParameter<double>("invisibleFraction")),
0027         exclusiveFraction_(conf.getParameter<double>("exclusiveFraction")),
0028         maxDistanceFilter_(conf.getParameter<bool>("maxDistanceFilter")),
0029         maxDistance_(conf.getParameter<double>("maxDistance")),
0030         maxDforTimingSquared_(conf.getParameter<double>("maxDforTimingSquared")),
0031         timeOffset_(conf.getParameter<double>("timeOffset")),
0032         minNHitsforTiming_(conf.getParameter<unsigned int>("minNHitsforTiming")),
0033         useMCFractionsForExclEnergy_(conf.getParameter<bool>("useMCFractionsForExclEnergy")),
0034         calibMinEta_(conf.getParameter<double>("calibMinEta")),
0035         calibMaxEta_(conf.getParameter<double>("calibMaxEta")),
0036         hadronCalib_(conf.getParameter<std::vector<double> >("hadronCalib")),
0037         egammaCalib_(conf.getParameter<std::vector<double> >("egammaCalib")),
0038         simClusterToken_(cc.consumes<SimClusterCollection>(conf.getParameter<edm::InputTag>("simClusterSrc"))),
0039         geomToken_(cc.esConsumes<edm::Transition::BeginRun>()) {}
0040
0041   ~RealisticSimClusterMapper() override {}
0042   RealisticSimClusterMapper(const RealisticSimClusterMapper&) = delete;
0043   RealisticSimClusterMapper& operator=(const RealisticSimClusterMapper&) = delete;
0044
0045   void updateEvent(const edm::Event&) final;
0046   void update(const edm::EventSetup&) final;
0047
0048   void buildClusters(const edm::Handle<reco::PFRecHitCollection>&,
0049                      const std::vector<bool>&,
0050                      const std::vector<bool>&,
0051                      reco::PFClusterCollection&) override;
0052
0053 private:
0054   hgcal::RecHitTools rhtools_;
0055   const float invisibleFraction_ = 0.3f;
0056   const float exclusiveFraction_ = 0.7f;
0057   const bool maxDistanceFilter_ = false;
0058   const float maxDistance_ = 10.f;
0059   const float maxDforTimingSquared_ = 4.0f;
0060   const float timeOffset_;
0061   const unsigned int minNHitsforTiming_ = 3;
0062   const bool useMCFractionsForExclEnergy_ = false;
0063   const float calibMinEta_ = 1.4;
0064   const float calibMaxEta_ = 3.0;
0065   std::vector<double> hadronCalib_;
0066   std::vector<double> egammaCalib_;
0067
0068   edm::EDGetTokenT<SimClusterCollection> simClusterToken_;
0069   edm::Handle<SimClusterCollection> simClusterH_;
0070
0071   edm::ESGetToken<CaloGeometry, CaloGeometryRecord> geomToken_;
0072 };
0073
0074 DEFINE_EDM_PLUGIN(InitialClusteringStepFactory, RealisticSimClusterMapper, "RealisticSimClusterMapper");
0075
0076 #ifdef PFLOW_DEBUG
0077 #define LOGVERB(x) edm::LogVerbatim(x)
0078 #define LOGWARN(x) edm::LogWarning(x)
0079 #define LOGERR(x) edm::LogError(x)
0080 #define LOGDRESSED(x) edm::LogInfo(x)
0081 #else
0082 #define LOGVERB(x) LogTrace(x)
0083 #define LOGWARN(x) edm::LogWarning(x)
0084 #define LOGERR(x) edm::LogError(x)
0085 #define LOGDRESSED(x) LogDebug(x)
0086 #endif
0087
0088 namespace {
0089
0090   inline bool isPi0(int pdgId) { return pdgId == 111; }
0091
0092   inline bool isEGamma(int pdgId) {
0093     pdgId = std::abs(pdgId);
0094     return (pdgId == 11) or (pdgId == 22);
0095   }
0096
0097   inline bool isHadron(int pdgId) {
0098     pdgId = std::abs(pdgId) % 10000;
0099     return ((pdgId > 100 and pdgId < 900) or (pdgId > 1000 and pdgId < 9000));
0100   }
0101 }  // namespace
0102
0103 void RealisticSimClusterMapper::updateEvent(const edm::Event& ev) { ev.getByToken(simClusterToken_, simClusterH_); }
0104
0105 void RealisticSimClusterMapper::update(const edm::EventSetup& es) { rhtools_.setGeometry(es.getData(geomToken_)); }
0106
0107 void RealisticSimClusterMapper::buildClusters(const edm::Handle<reco::PFRecHitCollection>& input,
0108                                               const std::vector<bool>& rechitMask,
0109                                               const std::vector<bool>& seedable,
0110                                               reco::PFClusterCollection& output) {
0111   const SimClusterCollection& simClusters = *simClusterH_;
0112   auto const& hits = *input;
0113   RealisticHitToClusterAssociator realisticAssociator;
0114   const int numberOfLayers = rhtools_.getGeometryType() == 0 ? rhtools_.getLayer(ForwardSubdetector::ForwardEmpty)
0115                                                              : rhtools_.getLayer(DetId::Forward);
0116   realisticAssociator.init(hits.size(), simClusters.size(), numberOfLayers + 1);
0117   // for quick indexing back to hit energy
0118   std::unordered_map<uint32_t, size_t> detIdToIndex(hits.size());
0119   for (uint32_t i = 0; i < hits.size(); ++i) {
0120     detIdToIndex[hits[i].detId()] = i;
0121     auto ref = makeRefhit(input, i);
0122     const auto& hitPos = rhtools_.getPosition(ref->detId());
0123
0124     realisticAssociator.insertHitPosition(hitPos.x(), hitPos.y(), hitPos.z(), i);
0125     realisticAssociator.insertHitEnergy(ref->energy(), i);
0126     realisticAssociator.insertLayerId(rhtools_.getLayerWithOffset(ref->detId()), i);
0127   }
0128
0129   for (unsigned int ic = 0; ic < simClusters.size(); ++ic) {
0130     const auto& sc = simClusters[ic];
0131     const auto& hitsAndFractions = sc.hits_and_fractions();
0132     for (const auto& hAndF : hitsAndFractions) {
0133       auto itr = detIdToIndex.find(hAndF.first);
0134       if (itr == detIdToIndex.end()) {
0135         continue;  // hit wasn't saved in reco or did not pass the SNR threshold
0136       }
0137       auto hitId = itr->second;
0138       auto ref = makeRefhit(input, hitId);
0139       float fraction = hAndF.second;
0140       float associatedEnergy = fraction * ref->energy();
0141       realisticAssociator.insertSimClusterIdAndFraction(ic, fraction, hitId, associatedEnergy);
0142     }
0143   }
0144   realisticAssociator.computeAssociation(
0145       exclusiveFraction_, useMCFractionsForExclEnergy_, rhtools_.lastLayerEE(), rhtools_.lastLayerFH());
0146   realisticAssociator.findAndMergeInvisibleClusters(invisibleFraction_, exclusiveFraction_);
0147   realisticAssociator.findCentersOfGravity();
0148   if (maxDistanceFilter_)
0149     realisticAssociator.filterHitsByDistance(maxDistance_);
0150
0151   const auto& realisticClusters = realisticAssociator.realisticClusters();
0152   unsigned int nClusters = realisticClusters.size();
0153   float bin_norm = 1. / (calibMaxEta_ - calibMinEta_);
0154   float egamma_bin_norm = egammaCalib_.size() * bin_norm;
0155   float hadron_bin_norm = hadronCalib_.size() * bin_norm;
0156   for (unsigned ic = 0; ic < nClusters; ++ic) {
0157     float highest_energy = 0.0f;
0158     output.emplace_back();
0159     reco::PFCluster& back = output.back();
0160     edm::Ref<std::vector<reco::PFRecHit> > seed;
0161     float energyCorrection = 1.f;
0162     float timeRealisticSC = -99.;
0163     if (realisticClusters[ic].isVisible()) {
0164       int pdgId = simClusters[ic].pdgId();
0165       auto abseta = std::abs(simClusters[ic].eta());
0166       if ((abseta >= calibMinEta_) and (abseta <= calibMaxEta_))  //protecting range
0167       {
0168         if ((isEGamma(pdgId) or isPi0(pdgId)) and !egammaCalib_.empty()) {
0169           unsigned int etabin = std::floor(((abseta - calibMinEta_) * egamma_bin_norm));
0170           energyCorrection = egammaCalib_[etabin];
0171         } else if (isHadron(pdgId) and !(isPi0(pdgId)) and
0172                    !hadronCalib_
0173                         .empty())  // this function is expensive.. should we treat as hadron everything which is not egamma?
0174         {
0175           unsigned int etabin = std::floor(((abseta - calibMinEta_) * hadron_bin_norm));
0176           energyCorrection = hadronCalib_[etabin];
0177         }
0178       }
0179       std::vector<float> timeHits;
0180       const auto& hitsIdsAndFractions = realisticClusters[ic].hitsIdsAndFractions();
0181       for (const auto& idAndF : hitsIdsAndFractions) {
0182         auto fraction = idAndF.second;
0183         if (fraction > 0.f) {
0184           auto ref = makeRefhit(input, idAndF.first);
0185           back.addRecHitFraction(reco::PFRecHitFraction(ref, fraction));
0186           const float hit_energy = fraction * ref->energy();
0187           if (hit_energy > highest_energy || highest_energy == 0.0) {
0188             highest_energy = hit_energy;
0189             seed = ref;
0190           }
0191           //select hits good for timing
0192           if (ref->time() > -1.) {
0193             int rhLayer = rhtools_.getLayerWithOffset(ref->detId());
0194             std::array<float, 3> scPosition = realisticClusters[ic].getCenterOfGravity(rhLayer);
0195             float distanceSquared =
0196                 std::pow((ref->position().x() - scPosition[0]), 2) + std::pow((ref->position().y() - scPosition[1]), 2);
0197             if (distanceSquared < maxDforTimingSquared_) {
0198               timeHits.push_back(ref->time());
0199             }
0200           }
0201         }
0202       }
0203       //assign time if minimum number of hits
0204       hgcalsimclustertime::ComputeClusterTime timeEstimator;
0205       timeRealisticSC = (timeEstimator.fixSizeHighestDensity(timeHits)).first;
0206     }
0207     if (!back.hitsAndFractions().empty()) {
0208       back.setSeed(seed->detId());
0209       back.setEnergy(realisticClusters[ic].getEnergy());
0210       back.setCorrectedEnergy(energyCorrection * realisticClusters[ic].getEnergy());  //applying energy correction
0211       back.setTime(timeRealisticSC);
0212     } else {
0213       back.setSeed(-1);
0214       back.setEnergy(0.f);
0215     }
0216   }
0217 }