CMS 3D CMS Logo

List of all members | Public Member Functions | Private Member Functions | Private Attributes
Basic2DGenericPFlowPositionCalc Class Reference
Inheritance diagram for Basic2DGenericPFlowPositionCalc:
PFCPositionCalculatorBase

Public Member Functions

 Basic2DGenericPFlowPositionCalc (const edm::ParameterSet &conf, edm::ConsumesCollector &cc)
 
 Basic2DGenericPFlowPositionCalc (const Basic2DGenericPFlowPositionCalc &)=delete
 
void calculateAndSetPosition (reco::PFCluster &, const HcalPFCuts *) override
 
void calculateAndSetPositions (reco::PFClusterCollection &, const HcalPFCuts *) override
 
Basic2DGenericPFlowPositionCalcoperator= (const Basic2DGenericPFlowPositionCalc &)=delete
 
- Public Member Functions inherited from PFCPositionCalculatorBase
const std::string & name () const
 
PosCalcoperator= (const PosCalc &)=delete
 
 PFCPositionCalculatorBase (const edm::ParameterSet &conf, edm::ConsumesCollector &cc)
 
 PFCPositionCalculatorBase (const PosCalc &)=delete
 
virtual void update (const edm::EventSetup &)
 
virtual ~PFCPositionCalculatorBase ()=default
 

Private Member Functions

void calculateAndSetPositionActual (reco::PFCluster &, const HcalPFCuts *) const
 

Private Attributes

std::tuple< std::vector< int >, std::vector< int >, std::vector< float > > _logWeightDenom
 
const float _minAllowedNorm
 
const int _posCalcNCrystals
 
std::unique_ptr< CaloRecHitResolutionProvider_timeResolutionCalcBarrel
 
std::unique_ptr< CaloRecHitResolutionProvider_timeResolutionCalcEndcap
 

Additional Inherited Members

- Protected Attributes inherited from PFCPositionCalculatorBase
const float _minFractionInCalc
 

Detailed Description

Definition at line 21 of file Basic2DGenericPFlowPositionCalc.cc.

Constructor & Destructor Documentation

◆ Basic2DGenericPFlowPositionCalc() [1/2]

Basic2DGenericPFlowPositionCalc::Basic2DGenericPFlowPositionCalc ( const edm::ParameterSet conf,
edm::ConsumesCollector cc 
)
inline

Definition at line 23 of file Basic2DGenericPFlowPositionCalc.cc.

References _logWeightDenom, _posCalcNCrystals, _timeResolutionCalcBarrel, _timeResolutionCalcEndcap, cms::cuda::assert(), HLT_2024v13_cff::depths, HLT_2024v13_cff::detectorEnum, Exception, edm::ParameterSet::getParameter(), edm::ParameterSet::getParameterSet(), edm::ParameterSet::getParameterSetVector(), mps_fire::i, muonDTDigis_cfi::pset, and AlCaHLTBitMon_QueryRunRegistry::string.

25  _posCalcNCrystals(conf.getParameter<int>("posCalcNCrystals")),
26  _minAllowedNorm(conf.getParameter<double>("minAllowedNormalization")) {
27  std::vector<int> detectorEnum;
28  std::vector<int> depths;
29  std::vector<double> logWeightDenom;
30  std::vector<float> logWeightDenomInv;
31 
32  const auto& logWeightDenominatorByDetectorPSet = conf.getParameterSetVector("logWeightDenominatorByDetector");
33  if (!logWeightDenominatorByDetectorPSet.empty()) {
34  for (const auto& pset : logWeightDenominatorByDetectorPSet) {
35  const auto& det = pset.getParameter<std::string>("detector");
36  if (det.empty()) {
37  throw cms::Exception("logWeightDenominatorByDetectorPSet") << "logWeightDenominator : detector not specified";
38  }
39 
40  if (det == std::string("HCAL_BARREL1") || det == std::string("HCAL_ENDCAP")) {
41  std::vector<int> depthsT = pset.getParameter<std::vector<int> >("depths");
42  std::vector<double> logWeightDenomT = pset.getParameter<std::vector<double> >("logWeightDenominator");
43  if (logWeightDenomT.size() != depthsT.size()) {
44  throw cms::Exception("logWeightDenominator") << "logWeightDenominator mismatch with the numbers of depths";
45  }
46  for (unsigned int i = 0; i < depthsT.size(); ++i) {
47  if (det == std::string("HCAL_BARREL1"))
48  detectorEnum.push_back(1);
49  if (det == std::string("HCAL_ENDCAP"))
50  detectorEnum.push_back(2);
51  depths.push_back(depthsT[i]);
52  logWeightDenom.push_back(logWeightDenomT[i]);
53  }
54  }
55  }
56  } else {
57  detectorEnum.push_back(0);
58  depths.push_back(0);
59  logWeightDenom.push_back(conf.getParameter<double>("logWeightDenominator"));
60  }
61 
62  for (unsigned int i = 0; i < depths.size(); ++i) {
63  logWeightDenomInv.push_back(1. / logWeightDenom[i]);
64  }
65 
66  // _logWeightDenom = std::make_pair(depths,logWeightDenomInv);
67  _logWeightDenom = std::make_tuple(detectorEnum, depths, logWeightDenomInv);
68 
69  _timeResolutionCalcBarrel.reset(nullptr);
70  const auto& timeResConfBarrel = conf.getParameterSet("timeResolutionCalcBarrel");
71  if (!timeResConfBarrel.empty() && timeResConfBarrel.getParameter<double>("threshHighE") >= 0)
72  _timeResolutionCalcBarrel = std::make_unique<CaloRecHitResolutionProvider>(timeResConfBarrel);
73  _timeResolutionCalcEndcap.reset(nullptr);
74  const auto& timeResConfEndcap = conf.getParameterSet("timeResolutionCalcEndcap");
75  if (!timeResConfEndcap.empty() && timeResConfEndcap.getParameter<double>("threshHighE") >= 0)
76  _timeResolutionCalcEndcap = std::make_unique<CaloRecHitResolutionProvider>(timeResConfEndcap);
77 
78  switch (_posCalcNCrystals) {
79  case 5:
80  case 9:
81  case -1:
82  break;
83  default:
84  edm::LogError("Basic2DGenericPFlowPositionCalc") << "posCalcNCrystals not valid";
85  assert(0); // bug
86  }
87  }
PFCPositionCalculatorBase(const edm::ParameterSet &conf, edm::ConsumesCollector &cc)
T getParameter(std::string const &) const
Definition: ParameterSet.h:307
uint32_t cc[maxCellsPerHit]
Definition: gpuFishbone.h:49
ParameterSet const & getParameterSet(std::string const &) const
Log< level::Error, false > LogError
assert(be >=bs)
std::unique_ptr< CaloRecHitResolutionProvider > _timeResolutionCalcBarrel
std::tuple< std::vector< int >, std::vector< int >, std::vector< float > > _logWeightDenom
std::unique_ptr< CaloRecHitResolutionProvider > _timeResolutionCalcEndcap
VParameterSet const & getParameterSetVector(std::string const &name) const

◆ Basic2DGenericPFlowPositionCalc() [2/2]

Basic2DGenericPFlowPositionCalc::Basic2DGenericPFlowPositionCalc ( const Basic2DGenericPFlowPositionCalc )
delete

Member Function Documentation

◆ calculateAndSetPosition()

void Basic2DGenericPFlowPositionCalc::calculateAndSetPosition ( reco::PFCluster cluster,
const HcalPFCuts hcalCuts 
)
overridevirtual

Implements PFCPositionCalculatorBase.

Definition at line 115 of file Basic2DGenericPFlowPositionCalc.cc.

References calculateAndSetPositionActual().

115  {
116  calculateAndSetPositionActual(cluster, hcalCuts);
117 }
void calculateAndSetPositionActual(reco::PFCluster &, const HcalPFCuts *) const

◆ calculateAndSetPositionActual()

void Basic2DGenericPFlowPositionCalc::calculateAndSetPositionActual ( reco::PFCluster cluster,
const HcalPFCuts hcalCuts 
) const
private

Definition at line 126 of file Basic2DGenericPFlowPositionCalc.cc.

References _logWeightDenom, _minAllowedNorm, PFCPositionCalculatorBase::_minFractionInCalc, _posCalcNCrystals, _timeResolutionCalcBarrel, _timeResolutionCalcEndcap, a, b, nano_mu_local_reco_cff::bool, reco::PFCluster::calculatePositionREP(), bookConverter::compute(), declareDynArray, hcalRecHitTable_cff::depth, reco::PFRecHit::depth(), HLT_2024v13_cff::detectorEnum, reco::PFRecHit::detId(), hcalRecHitTable_cff::energy, Exception, f, myMath::fast_logf(), nano_mu_digi_cff::float, HLT_2024v13_cff::fraction, HcalCondObjectContainer< Item >::getValues(), PFLayer::HCAL_BARREL1, PFLayer::HCAL_ENDCAP, photonIsolationHIProducer_cfi::hf, hfClusterShapes_cfi::hits, mps_fire::i, createfilelist::int, PixelPluginsPhase0_cfi::isBarrel, edm::isNotFinite(), B2GTnPMonitor_cfi::item, dqmiolumiharvest::j, dqmdumpme::k, reco::PFRecHit::layer(), SiStripPI::max, ecalTB2006H4_GenSimDigiReco_cfg::mySeed, TrackingDataMCValidation_Standalone_cff::nhits, PFLayer::NONE, AlCaHLTBitMon_ParallelJobs::p, reco::PFRecHit::position(), DetId::rawId(), reco::PFCluster::recHitFractions(), reco::CaloCluster::seed(), DBoxMetadataHelper::set_intersection(), reco::PFCluster::setDepth(), reco::CaloCluster::setEnergy(), reco::PFCluster::setLayer(), reco::CaloCluster::setPosition(), reco::PFCluster::setTime(), reco::PFCluster::setTimeError(), jetUpdater_cfi::sort, mathSSE::sqrt(), DiMuonV_cfg::threshold, reco::PFRecHit::time(), mitigatedMETSequence_cff::U, unInitDynArray, UNLIKELY, x, y, and z.

Referenced by calculateAndSetPosition(), and calculateAndSetPositions().

127  {
128  if (!cluster.seed()) {
129  throw cms::Exception("ClusterWithNoSeed") << " Found a cluster with no seed: " << cluster;
130  }
131  double cl_energy = 0;
132  double cl_time = 0;
133  double cl_timeweight = 0.0;
134  double max_e = 0.0;
135  PFLayer::Layer max_e_layer = PFLayer::NONE;
136  // find the seed and max layer and also calculate time
137  //Michalis : Even if we dont use timing in clustering here we should fill
138  //the time information for the cluster. This should use the timing resolution(1/E)
139  //so the weight should be fraction*E^2
140  //calculate a simplistic depth now. The log weighted will be done
141  //in different stage
142 
143  auto const recHitCollection =
144  &(*cluster.recHitFractions()[0].recHitRef()) - cluster.recHitFractions()[0].recHitRef().key();
145  auto nhits = cluster.recHitFractions().size();
146  struct LHit {
147  reco::PFRecHit const* hit;
148  float energy;
149  float fraction;
150  };
151  declareDynArray(LHit, nhits, hits);
152  for (auto i = 0U; i < nhits; ++i) {
153  auto const& hf = cluster.recHitFractions()[i];
154  auto k = hf.recHitRef().key();
155  auto p = recHitCollection + k;
156  hits[i] = {p, (*p).energy(), float(hf.fraction())};
157  }
158 
160  LHit mySeed = {};
161  for (auto const& rhf : hits) {
162  const reco::PFRecHit& refhit = *rhf.hit;
163  if (refhit.detId() == cluster.seed())
164  mySeed = rhf;
165  const auto rh_fraction = rhf.fraction;
166  const auto rh_rawenergy = rhf.energy;
167  const auto rh_energy = rh_rawenergy * rh_fraction;
168 #ifdef PF_DEBUG
169  if UNLIKELY (edm::isNotFinite(rh_energy)) {
170  throw cms::Exception("PFClusterAlgo") << "rechit " << refhit.detId() << " has a NaN energy... "
171  << "The input of the particle flow clustering seems to be corrupted.";
172  }
173 #endif
174  cl_energy += rh_energy;
175  // If time resolution is given, calculated weighted average
176  if (resGiven) {
177  double res2 = 1.e-4;
178  int cell_layer = (int)refhit.layer();
179  res2 = isBarrel(cell_layer) ? 1. / _timeResolutionCalcBarrel->timeResolution2(rh_rawenergy)
180  : 1. / _timeResolutionCalcEndcap->timeResolution2(rh_rawenergy);
181  cl_time += rh_fraction * refhit.time() * res2;
182  cl_timeweight += rh_fraction * res2;
183  } else { // assume resolution = 1/E**2
184  const double rh_rawenergy2 = rh_rawenergy * rh_rawenergy;
185  cl_timeweight += rh_rawenergy2 * rh_fraction;
186  cl_time += rh_rawenergy2 * rh_fraction * refhit.time();
187  }
188 
189  if (rh_energy > max_e) {
190  max_e = rh_energy;
191  max_e_layer = refhit.layer();
192  }
193  }
194 
195  cluster.setEnergy(cl_energy);
196  cluster.setTime(cl_time / cl_timeweight);
197  if (resGiven) {
198  cluster.setTimeError(std::sqrt(1.0f / float(cl_timeweight)));
199  }
200  cluster.setLayer(max_e_layer);
201 
202  // calculate the position
203  bool single_depth = true;
204  int ref_depth = -1;
205  double depth = 0.0;
206  double position_norm = 0.0;
207  double x(0.0), y(0.0), z(0.0);
208  if (nullptr != mySeed.hit) {
209  auto seedNeighbours = mySeed.hit->neighbours();
210  switch (_posCalcNCrystals) {
211  case 5:
212  seedNeighbours = mySeed.hit->neighbours4();
213  break;
214  case 9:
215  seedNeighbours = mySeed.hit->neighbours8();
216  break;
217  default:
218  break;
219  }
220 
221  auto compute = [&](LHit const& rhf) {
222  const reco::PFRecHit& refhit = *rhf.hit;
223 
224  int cell_layer = (int)refhit.layer();
225  float threshold = 0;
226 
227  if (hcalCuts != nullptr && // this means, cutsFromDB is set to True in the producer code
228  (cell_layer == PFLayer::HCAL_BARREL1 || cell_layer == PFLayer::HCAL_ENDCAP)) {
229  HcalDetId thisId = refhit.detId();
230  const HcalPFCut* item = hcalCuts->getValues(thisId.rawId());
231  threshold = 1. / (item->noiseThreshold());
232 
233  } else {
234  for (unsigned int j = 0; j < (std::get<2>(_logWeightDenom)).size(); ++j) {
235  // barrel is detecor type1
236  int detectorEnum = std::get<0>(_logWeightDenom)[j];
237  int depth = std::get<1>(_logWeightDenom)[j];
238  if ((cell_layer == PFLayer::HCAL_BARREL1 && detectorEnum == 1 && refhit.depth() == depth) ||
239  (cell_layer == PFLayer::HCAL_ENDCAP && detectorEnum == 2 && refhit.depth() == depth) ||
240  detectorEnum == 0) {
241  threshold = std::get<2>(_logWeightDenom)[j];
242  }
243  }
244  }
245 
246  if (ref_depth < 0)
247  ref_depth = refhit.depth(); // Initialize reference depth
248  else if (refhit.depth() != ref_depth) {
249  // Found a rechit with a different depth
250  single_depth = false;
251  }
252  const auto rh_energy = rhf.energy * rhf.fraction;
253  const auto norm =
254  (rhf.fraction < _minFractionInCalc ? 0.0f : std::max(0.0f, vdt::fast_logf(rh_energy * threshold)));
255  const auto rhpos_xyz = refhit.position() * norm;
256  x += rhpos_xyz.x();
257  y += rhpos_xyz.y();
258  z += rhpos_xyz.z();
259  depth += refhit.depth() * norm;
260  position_norm += norm;
261  };
262 
263  if (_posCalcNCrystals != -1) // sorted to make neighbour search faster (maybe)
264  std::sort(hits.begin(), hits.end(), [](LHit const& a, LHit const& b) { return a.hit < b.hit; });
265 
266  if (_posCalcNCrystals == -1)
267  for (auto const& rhf : hits)
268  compute(rhf);
269  else { // only seed and its neighbours
270  compute(mySeed);
271  // search seedNeighbours to find energy fraction in cluster (sic)
272  unInitDynArray(reco::PFRecHit const*, seedNeighbours.size(), nei);
273  for (auto k : seedNeighbours) {
274  nei.push_back(recHitCollection + k);
275  }
276  std::sort(nei.begin(), nei.end());
277  struct LHitLess {
278  auto operator()(LHit const& a, reco::PFRecHit const* b) const { return a.hit < b; }
279  auto operator()(reco::PFRecHit const* b, LHit const& a) const { return b < a.hit; }
280  };
282  hits.begin(), hits.end(), nei.begin(), nei.end(), boost::make_function_output_iterator(compute), LHitLess());
283  }
284  } else {
285  throw cms::Exception("Basic2DGenerticPFlowPositionCalc")
286  << "Cluster seed hit is null, something is wrong with PFlow RecHit!";
287  }
288 
289  if (position_norm < _minAllowedNorm) {
290  edm::LogError("WeirdClusterNormalization") << "PFCluster too far from seeding cell: set position to (0,0,0).";
291  cluster.setPosition(math::XYZPoint(0, 0, 0));
292  cluster.calculatePositionREP();
293  } else {
294  const double norm_inverse = 1.0 / position_norm;
295  x *= norm_inverse;
296  y *= norm_inverse;
297  z *= norm_inverse;
298  if (single_depth)
299  depth = ref_depth;
300  else
301  depth *= norm_inverse;
302  cluster.setPosition(math::XYZPoint(x, y, z));
303  cluster.setDepth(depth);
304  cluster.calculatePositionREP();
305  }
306 }
void setLayer(PFLayer::Layer layer)
set layer
Definition: PFCluster.cc:49
constexpr bool isNotFinite(T x)
Definition: isFinite.h:9
void setPosition(const math::XYZPoint &p)
Definition: CaloCluster.h:140
DetId seed() const
return DetId of seed
Definition: CaloCluster.h:219
void setTime(float time, float timeError=0)
Definition: PFCluster.h:84
PFLayer::Layer layer() const
rechit layer
Definition: PFRecHit.h:96
const std::vector< reco::PFRecHitFraction > & recHitFractions() const
vector of rechit fractions
Definition: PFCluster.h:65
void setEnergy(double energy)
Definition: CaloCluster.h:136
float time() const
timing for cleaned hits
Definition: PFRecHit.h:102
Log< level::Error, false > LogError
PositionType const & position() const
rechit cell centre x, y, z
Definition: PFRecHit.h:117
const Item * getValues(DetId fId, bool throwOnFail=true) const
unsigned detId() const
rechit detId
Definition: PFRecHit.h:93
Particle flow rechit (rechit + geometry and topology information). See clustering algorithm in PFClus...
Definition: PFRecHit.h:31
void calculatePositionREP()
computes posrep_ once and for all
Definition: PFCluster.h:95
T sqrt(T t)
Definition: SSEVec.h:23
#define unInitDynArray(T, n, x)
Definition: DynArray.h:88
double f[11][100]
Layer
layer definition
Definition: PFLayer.h:29
void setTimeError(float timeError)
Definition: PFCluster.h:88
constexpr uint32_t rawId() const
get the raw id
Definition: DetId.h:57
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
def compute(min, max)
std::unique_ptr< CaloRecHitResolutionProvider > _timeResolutionCalcBarrel
double b
Definition: hdecay.h:120
std::tuple< std::vector< int >, std::vector< int >, std::vector< float > > _logWeightDenom
std::unique_ptr< CaloRecHitResolutionProvider > _timeResolutionCalcEndcap
double a
Definition: hdecay.h:121
void setDepth(double depth)
Definition: PFCluster.h:89
#define declareDynArray(T, n, x)
Definition: DynArray.h:91
int depth() const
depth for segemntation
Definition: PFRecHit.h:105
#define UNLIKELY(x)
Definition: Likely.h:21
float fast_logf(float x)
std::vector< std::string > set_intersection(std::vector< std::string > const &v1, std::vector< std::string > const &v2)

◆ calculateAndSetPositions()

void Basic2DGenericPFlowPositionCalc::calculateAndSetPositions ( reco::PFClusterCollection clusters,
const HcalPFCuts hcalCuts 
)
overridevirtual

Implements PFCPositionCalculatorBase.

Definition at line 119 of file Basic2DGenericPFlowPositionCalc.cc.

References calculateAndSetPositionActual(), and bsc_activity_cfg::clusters.

120  {
121  for (reco::PFCluster& cluster : clusters) {
122  calculateAndSetPositionActual(cluster, hcalCuts);
123  }
124 }
Particle flow cluster, see clustering algorithm in PFClusterAlgo.
Definition: PFCluster.h:42
void calculateAndSetPositionActual(reco::PFCluster &, const HcalPFCuts *) const

◆ operator=()

Basic2DGenericPFlowPositionCalc& Basic2DGenericPFlowPositionCalc::operator= ( const Basic2DGenericPFlowPositionCalc )
delete

Member Data Documentation

◆ _logWeightDenom

std::tuple<std::vector<int>, std::vector<int>, std::vector<float> > Basic2DGenericPFlowPositionCalc::_logWeightDenom
private

◆ _minAllowedNorm

const float Basic2DGenericPFlowPositionCalc::_minAllowedNorm
private

Definition at line 98 of file Basic2DGenericPFlowPositionCalc.cc.

Referenced by calculateAndSetPositionActual().

◆ _posCalcNCrystals

const int Basic2DGenericPFlowPositionCalc::_posCalcNCrystals
private

◆ _timeResolutionCalcBarrel

std::unique_ptr<CaloRecHitResolutionProvider> Basic2DGenericPFlowPositionCalc::_timeResolutionCalcBarrel
private

◆ _timeResolutionCalcEndcap

std::unique_ptr<CaloRecHitResolutionProvider> Basic2DGenericPFlowPositionCalc::_timeResolutionCalcEndcap
private