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_2024v11_cff::depths, HLT_2024v11_cff::detectorEnum, Exception, edm::ParameterSet::exists(), 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  if (conf.exists("logWeightDenominatorByDetector")) {
33  const std::vector<edm::ParameterSet>& logWeightDenominatorByDetectorPSet =
34  conf.getParameterSetVector("logWeightDenominatorByDetector");
35 
36  for (const auto& pset : logWeightDenominatorByDetectorPSet) {
37  if (!pset.exists("detector")) {
38  throw cms::Exception("logWeightDenominatorByDetectorPSet") << "logWeightDenominator : detector not specified";
39  }
40 
41  const std::string& det = pset.getParameter<std::string>("detector");
42 
43  if (det == std::string("HCAL_BARREL1") || det == std::string("HCAL_ENDCAP")) {
44  std::vector<int> depthsT = pset.getParameter<std::vector<int> >("depths");
45  std::vector<double> logWeightDenomT = pset.getParameter<std::vector<double> >("logWeightDenominator");
46  if (logWeightDenomT.size() != depthsT.size()) {
47  throw cms::Exception("logWeightDenominator") << "logWeightDenominator mismatch with the numbers of depths";
48  }
49  for (unsigned int i = 0; i < depthsT.size(); ++i) {
50  if (det == std::string("HCAL_BARREL1"))
51  detectorEnum.push_back(1);
52  if (det == std::string("HCAL_ENDCAP"))
53  detectorEnum.push_back(2);
54  depths.push_back(depthsT[i]);
55  logWeightDenom.push_back(logWeightDenomT[i]);
56  }
57  }
58  }
59  } else {
60  detectorEnum.push_back(0);
61  depths.push_back(0);
62  logWeightDenom.push_back(conf.getParameter<double>("logWeightDenominator"));
63  }
64 
65  for (unsigned int i = 0; i < depths.size(); ++i) {
66  logWeightDenomInv.push_back(1. / logWeightDenom[i]);
67  }
68 
69  // _logWeightDenom = std::make_pair(depths,logWeightDenomInv);
70  _logWeightDenom = std::make_tuple(detectorEnum, depths, logWeightDenomInv);
71 
72  _timeResolutionCalcBarrel.reset(nullptr);
73  if (conf.exists("timeResolutionCalcBarrel")) {
74  const edm::ParameterSet& timeResConf = conf.getParameterSet("timeResolutionCalcBarrel");
75  _timeResolutionCalcBarrel = std::make_unique<CaloRecHitResolutionProvider>(timeResConf);
76  }
77  _timeResolutionCalcEndcap.reset(nullptr);
78  if (conf.exists("timeResolutionCalcEndcap")) {
79  const edm::ParameterSet& timeResConf = conf.getParameterSet("timeResolutionCalcEndcap");
80  _timeResolutionCalcEndcap = std::make_unique<CaloRecHitResolutionProvider>(timeResConf);
81  }
82 
83  switch (_posCalcNCrystals) {
84  case 5:
85  case 9:
86  case -1:
87  break;
88  default:
89  edm::LogError("Basic2DGenericPFlowPositionCalc") << "posCalcNCrystals not valid";
90  assert(0); // bug
91  }
92  }
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
bool exists(std::string const &parameterName) const
checks if a parameter exists
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 120 of file Basic2DGenericPFlowPositionCalc.cc.

References calculateAndSetPositionActual().

120  {
121  calculateAndSetPositionActual(cluster, hcalCuts);
122 }
void calculateAndSetPositionActual(reco::PFCluster &, const HcalPFCuts *) const

◆ calculateAndSetPositionActual()

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

Definition at line 131 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_2024v11_cff::detectorEnum, reco::PFRecHit::detId(), hcalRecHitTable_cff::energy, Exception, f, myMath::fast_logf(), nano_mu_digi_cff::float, HLT_2024v11_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, 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().

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

References calculateAndSetPositionActual(), and bsc_activity_cfg::clusters.

125  {
126  for (reco::PFCluster& cluster : clusters) {
127  calculateAndSetPositionActual(cluster, hcalCuts);
128  }
129 }
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 103 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