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