CMS 3D CMS Logo

RealisticSimClusterMapper.cc
Go to the documentation of this file.
1 // Author: Felice Pantaleo
3 // Date: 30/06/2017
4 // Email: felice@cern.ch
6 
7 #include "RealisticCluster.h"
8 
21 
22 #include <unordered_map>
23 
25 public:
28  invisibleFraction_(conf.getParameter<double>("invisibleFraction")),
29  exclusiveFraction_(conf.getParameter<double>("exclusiveFraction")),
30  maxDistanceFilter_(conf.getParameter<bool>("maxDistanceFilter")),
31  maxDistance_(conf.getParameter<double>("maxDistance")),
32  maxDforTimingSquared_(conf.getParameter<double>("maxDforTimingSquared")),
33  timeOffset_(conf.getParameter<double>("timeOffset")),
34  minNHitsforTiming_(conf.getParameter<unsigned int>("minNHitsforTiming")),
35  useMCFractionsForExclEnergy_(conf.getParameter<bool>("useMCFractionsForExclEnergy")),
36  calibMinEta_(conf.getParameter<double>("calibMinEta")),
37  calibMaxEta_(conf.getParameter<double>("calibMaxEta")),
38  hadronCalib_(conf.getParameter<std::vector<double> >("hadronCalib")),
39  egammaCalib_(conf.getParameter<std::vector<double> >("egammaCalib")),
40  simClusterToken_(cc.consumes<SimClusterCollection>(conf.getParameter<edm::InputTag>("simClusterSrc"))),
41  geomToken_(cc.esConsumes<edm::Transition::BeginRun>()) {}
42 
46 
47  void updateEvent(const edm::Event&) final;
48  void update(const edm::EventSetup&) final;
49 
51  const std::vector<bool>&,
52  const std::vector<bool>&,
54  const HcalPFCuts*) override;
55 
56 private:
58  const float invisibleFraction_ = 0.3f;
59  const float exclusiveFraction_ = 0.7f;
60  const bool maxDistanceFilter_ = false;
61  const float maxDistance_ = 10.f;
62  const float maxDforTimingSquared_ = 4.0f;
63  const float timeOffset_;
64  const unsigned int minNHitsforTiming_ = 3;
65  const bool useMCFractionsForExclEnergy_ = false;
66  const float calibMinEta_ = 1.4;
67  const float calibMaxEta_ = 3.0;
68  std::vector<double> hadronCalib_;
69  std::vector<double> egammaCalib_;
70 
73 
75 };
76 
78 
79 #ifdef PFLOW_DEBUG
80 #define LOGVERB(x) edm::LogVerbatim(x)
81 #define LOGWARN(x) edm::LogWarning(x)
82 #define LOGERR(x) edm::LogError(x)
83 #define LOGDRESSED(x) edm::LogInfo(x)
84 #else
85 #define LOGVERB(x) LogTrace(x)
86 #define LOGWARN(x) edm::LogWarning(x)
87 #define LOGERR(x) edm::LogError(x)
88 #define LOGDRESSED(x) LogDebug(x)
89 #endif
90 
91 namespace {
92 
93  inline bool isPi0(int pdgId) { return pdgId == 111; }
94 
95  inline bool isEGamma(int pdgId) {
96  pdgId = std::abs(pdgId);
97  return (pdgId == 11) or (pdgId == 22);
98  }
99 
100  inline bool isHadron(int pdgId) {
101  pdgId = std::abs(pdgId) % 10000;
102  return ((pdgId > 100 and pdgId < 900) or (pdgId > 1000 and pdgId < 9000));
103  }
104 } // namespace
105 
107 
109 
111  const std::vector<bool>& rechitMask,
112  const std::vector<bool>& seedable,
114  const HcalPFCuts* hcalCuts) {
115  const SimClusterCollection& simClusters = *simClusterH_;
116  auto const& hits = *input;
117  RealisticHitToClusterAssociator realisticAssociator;
120  realisticAssociator.init(hits.size(), simClusters.size(), numberOfLayers + 1);
121  // for quick indexing back to hit energy
122  std::unordered_map<uint32_t, size_t> detIdToIndex(hits.size());
123  for (uint32_t i = 0; i < hits.size(); ++i) {
124  detIdToIndex[hits[i].detId()] = i;
125  auto ref = makeRefhit(input, i);
126  const auto& hitPos = rhtools_.getPosition(ref->detId());
127 
128  realisticAssociator.insertHitPosition(hitPos.x(), hitPos.y(), hitPos.z(), i);
129  realisticAssociator.insertHitEnergy(ref->energy(), i);
130  realisticAssociator.insertLayerId(rhtools_.getLayerWithOffset(ref->detId()), i);
131  }
132 
133  for (unsigned int ic = 0; ic < simClusters.size(); ++ic) {
134  const auto& sc = simClusters[ic];
135  const auto& hitsAndFractions = sc.hits_and_fractions();
136  for (const auto& hAndF : hitsAndFractions) {
137  auto itr = detIdToIndex.find(hAndF.first);
138  if (itr == detIdToIndex.end()) {
139  continue; // hit wasn't saved in reco or did not pass the SNR threshold
140  }
141  auto hitId = itr->second;
142  auto ref = makeRefhit(input, hitId);
143  float fraction = hAndF.second;
144  float associatedEnergy = fraction * ref->energy();
145  realisticAssociator.insertSimClusterIdAndFraction(ic, fraction, hitId, associatedEnergy);
146  }
147  }
148  realisticAssociator.computeAssociation(
151  realisticAssociator.findCentersOfGravity();
152  if (maxDistanceFilter_)
153  realisticAssociator.filterHitsByDistance(maxDistance_);
154 
155  const auto& realisticClusters = realisticAssociator.realisticClusters();
156  unsigned int nClusters = realisticClusters.size();
157  float bin_norm = 1. / (calibMaxEta_ - calibMinEta_);
158  float egamma_bin_norm = egammaCalib_.size() * bin_norm;
159  float hadron_bin_norm = hadronCalib_.size() * bin_norm;
160  for (unsigned ic = 0; ic < nClusters; ++ic) {
161  float highest_energy = 0.0f;
162  output.emplace_back();
163  reco::PFCluster& back = output.back();
165  float energyCorrection = 1.f;
166  float timeRealisticSC = -99.;
167  if (realisticClusters[ic].isVisible()) {
168  int pdgId = simClusters[ic].pdgId();
169  auto abseta = std::abs(simClusters[ic].eta());
170  if ((abseta >= calibMinEta_) and (abseta <= calibMaxEta_)) //protecting range
171  {
172  if ((isEGamma(pdgId) or isPi0(pdgId)) and !egammaCalib_.empty()) {
173  unsigned int etabin = std::floor(((abseta - calibMinEta_) * egamma_bin_norm));
174  energyCorrection = egammaCalib_[etabin];
175  } else if (isHadron(pdgId) and !(isPi0(pdgId)) and
176  !hadronCalib_
177  .empty()) // this function is expensive.. should we treat as hadron everything which is not egamma?
178  {
179  unsigned int etabin = std::floor(((abseta - calibMinEta_) * hadron_bin_norm));
180  energyCorrection = hadronCalib_[etabin];
181  }
182  }
183  std::vector<float> timeHits;
184  const auto& hitsIdsAndFractions = realisticClusters[ic].hitsIdsAndFractions();
185  for (const auto& idAndF : hitsIdsAndFractions) {
186  auto fraction = idAndF.second;
187  if (fraction > 0.f) {
188  auto ref = makeRefhit(input, idAndF.first);
190  const float hit_energy = fraction * ref->energy();
191  if (hit_energy > highest_energy || highest_energy == 0.0) {
192  highest_energy = hit_energy;
193  seed = ref;
194  }
195  //select hits good for timing
196  if (ref->time() > -1.) {
197  int rhLayer = rhtools_.getLayerWithOffset(ref->detId());
198  std::array<float, 3> scPosition = realisticClusters[ic].getCenterOfGravity(rhLayer);
199  float distanceSquared =
200  std::pow((ref->position().x() - scPosition[0]), 2) + std::pow((ref->position().y() - scPosition[1]), 2);
201  if (distanceSquared < maxDforTimingSquared_) {
202  timeHits.push_back(ref->time());
203  }
204  }
205  }
206  }
207  //assign time if minimum number of hits
209  timeRealisticSC = (timeEstimator.fixSizeHighestDensity(timeHits)).first;
210  }
211  if (!back.hitsAndFractions().empty()) {
212  back.setSeed(seed->detId());
213  back.setEnergy(realisticClusters[ic].getEnergy());
214  back.setCorrectedEnergy(energyCorrection * realisticClusters[ic].getEnergy()); //applying energy correction
215  back.setTime(timeRealisticSC);
216  } else {
217  back.setSeed(-1);
218  back.setEnergy(0.f);
219  }
220  }
221 }
void insertLayerId(unsigned int layerId, unsigned int hitIndex)
ESGetTokenH3DDVariant esConsumes(std::string const &Record, edm::ConsumesCollector &)
Definition: DeDxTools.cc:283
const std::vector< std::pair< DetId, float > > & hitsAndFractions() const
Definition: CaloCluster.h:209
reco::PFRecHitRef makeRefhit(const edm::Handle< reco::PFRecHitCollection > &h, const unsigned i) const
T const & getData(const ESGetToken< T, R > &iToken) const noexcept(false)
Definition: EventSetup.h:119
unsigned int getLayer(DetId::Detector type, bool nose=false) const
Definition: RecHitTools.cc:318
void buildClusters(const edm::Handle< reco::PFRecHitCollection > &, const std::vector< bool > &, const std::vector< bool > &, reco::PFClusterCollection &, const HcalPFCuts *) override
Particle flow cluster, see clustering algorithm in PFClusterAlgo.
Definition: PFCluster.h:42
constexpr uint32_t numberOfLayers
uint32_t cc[maxCellsPerHit]
Definition: gpuFishbone.h:49
double getEnergy(HBHERecHitCollection::const_iterator hit, int useRaw=0, bool debug=false)
void setTime(float time, float timeError=0)
Definition: PFCluster.h:84
Fraction of a PFRecHit (rechits can be shared between several PFCluster&#39;s)
void setEnergy(double energy)
Definition: CaloCluster.h:135
int getGeometryType() const
Definition: RecHitTools.h:86
static std::string const input
Definition: EdmProvDump.cc:50
void setSeed(const DetId &id)
Definition: CaloCluster.h:145
RealisticSimClusterMapper(const edm::ParameterSet &conf, edm::ConsumesCollector &cc)
std::pair< float, float > fixSizeHighestDensity(std::vector< float > &time, std::vector< float > weight=std::vector< float >(), unsigned int minNhits=3, float deltaT=0.210, float timeWidthBy=0.5)
void init(std::size_t numberOfHits, std::size_t numberOfSimClusters, std::size_t numberOfLayers)
unsigned int lastLayerFH() const
Definition: RecHitTools.h:77
GlobalPoint getPosition(const DetId &id) const
Definition: RecHitTools.cc:140
const std::vector< RealisticCluster > & realisticClusters() const
void setCorrectedEnergy(double cenergy)
Definition: CaloCluster.h:136
void insertSimClusterIdAndFraction(unsigned int scIdx, float fraction, unsigned int hitIndex, float associatedEnergy)
void computeAssociation(float exclusiveFraction, bool useMCFractionsForExclEnergy, unsigned int fhOffset, unsigned int bhOffset)
The Signals That Services Can Subscribe To This is based on ActivityRegistry and is current per Services can connect to the signals distributed by the ActivityRegistry in order to monitor the activity of the application Each possible callback has some defined which we here list in angle e< void, edm::EventID const &, edm::Timestamp const & > We also list in braces which AR_WATCH_USING_METHOD_ is used for those or
Definition: Activities.doc:12
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
Transition
Definition: Transition.h:12
void findAndMergeInvisibleClusters(float invisibleFraction, float exclusiveFraction)
void insertHitPosition(float x, float y, float z, unsigned int hitIndex)
double f[11][100]
RealisticSimClusterMapper & operator=(const RealisticSimClusterMapper &)=delete
void update(const edm::EventSetup &) final
edm::EDGetTokenT< SimClusterCollection > simClusterToken_
void insertHitEnergy(float energy, unsigned int hitIndex)
void updateEvent(const edm::Event &) final
void addRecHitFraction(const reco::PFRecHitFraction &frac)
add a given fraction of the rechit
Definition: PFCluster.cc:33
void setGeometry(CaloGeometry const &)
Definition: RecHitTools.cc:79
HLT enums.
std::vector< PFCluster > PFClusterCollection
collection of PFCluster objects
Definition: PFClusterFwd.h:9
#define DEFINE_EDM_PLUGIN(factory, type, name)
Definition: output.py:1
edm::Handle< SimClusterCollection > simClusterH_
std::vector< SimCluster > SimClusterCollection
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
unsigned int lastLayerEE(bool nose=false) const
Definition: RecHitTools.h:76
edm::ESGetToken< CaloGeometry, CaloGeometryRecord > geomToken_
unsigned int getLayerWithOffset(const DetId &) const
Definition: RecHitTools.cc:376