File indexing completed on 2023-03-17 11:21:02
0001 #include "CommonTools/Utils/interface/DynArray.h"
0002 #include "DataFormats/ParticleFlowReco/interface/PFRecHit.h"
0003 #include "DataFormats/ParticleFlowReco/interface/PFRecHitFwd.h"
0004 #include "FWCore/Framework/interface/ConsumesCollector.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 #include "FWCore/Utilities/interface/isFinite.h"
0007 #include "RecoParticleFlow/PFClusterProducer/interface/CaloRecHitResolutionProvider.h"
0008 #include "RecoParticleFlow/PFClusterProducer/interface/PFCPositionCalculatorBase.h"
0009
0010 #include "vdt/vdtMath.h"
0011
0012 #include <boost/iterator/function_output_iterator.hpp>
0013
0014 #include <cmath>
0015 #include <iterator>
0016 #include <memory>
0017 #include <tuple>
0018
0019 class Basic2DGenericPFlowPositionCalc : public PFCPositionCalculatorBase {
0020 public:
0021 Basic2DGenericPFlowPositionCalc(const edm::ParameterSet& conf, edm::ConsumesCollector& cc)
0022 : PFCPositionCalculatorBase(conf, cc),
0023 _posCalcNCrystals(conf.getParameter<int>("posCalcNCrystals")),
0024 _minAllowedNorm(conf.getParameter<double>("minAllowedNormalization")) {
0025 std::vector<int> detectorEnum;
0026 std::vector<int> depths;
0027 std::vector<double> logWeightDenom;
0028 std::vector<float> logWeightDenomInv;
0029
0030 if (conf.exists("logWeightDenominatorByDetector")) {
0031 const std::vector<edm::ParameterSet>& logWeightDenominatorByDetectorPSet =
0032 conf.getParameterSetVector("logWeightDenominatorByDetector");
0033
0034 for (const auto& pset : logWeightDenominatorByDetectorPSet) {
0035 if (!pset.exists("detector")) {
0036 throw cms::Exception("logWeightDenominatorByDetectorPSet") << "logWeightDenominator : detector not specified";
0037 }
0038
0039 const std::string& det = pset.getParameter<std::string>("detector");
0040
0041 if (det == std::string("HCAL_BARREL1") || det == std::string("HCAL_ENDCAP")) {
0042 std::vector<int> depthsT = pset.getParameter<std::vector<int> >("depths");
0043 std::vector<double> logWeightDenomT = pset.getParameter<std::vector<double> >("logWeightDenominator");
0044 if (logWeightDenomT.size() != depthsT.size()) {
0045 throw cms::Exception("logWeightDenominator") << "logWeightDenominator mismatch with the numbers of depths";
0046 }
0047 for (unsigned int i = 0; i < depthsT.size(); ++i) {
0048 if (det == std::string("HCAL_BARREL1"))
0049 detectorEnum.push_back(1);
0050 if (det == std::string("HCAL_ENDCAP"))
0051 detectorEnum.push_back(2);
0052 depths.push_back(depthsT[i]);
0053 logWeightDenom.push_back(logWeightDenomT[i]);
0054 }
0055 }
0056 }
0057 } else {
0058 detectorEnum.push_back(0);
0059 depths.push_back(0);
0060 logWeightDenom.push_back(conf.getParameter<double>("logWeightDenominator"));
0061 }
0062
0063 for (unsigned int i = 0; i < depths.size(); ++i) {
0064 logWeightDenomInv.push_back(1. / logWeightDenom[i]);
0065 }
0066
0067
0068 _logWeightDenom = std::make_tuple(detectorEnum, depths, logWeightDenomInv);
0069
0070 _timeResolutionCalcBarrel.reset(nullptr);
0071 if (conf.exists("timeResolutionCalcBarrel")) {
0072 const edm::ParameterSet& timeResConf = conf.getParameterSet("timeResolutionCalcBarrel");
0073 _timeResolutionCalcBarrel = std::make_unique<CaloRecHitResolutionProvider>(timeResConf);
0074 }
0075 _timeResolutionCalcEndcap.reset(nullptr);
0076 if (conf.exists("timeResolutionCalcEndcap")) {
0077 const edm::ParameterSet& timeResConf = conf.getParameterSet("timeResolutionCalcEndcap");
0078 _timeResolutionCalcEndcap = std::make_unique<CaloRecHitResolutionProvider>(timeResConf);
0079 }
0080
0081 switch (_posCalcNCrystals) {
0082 case 5:
0083 case 9:
0084 case -1:
0085 break;
0086 default:
0087 edm::LogError("Basic2DGenericPFlowPositionCalc") << "posCalcNCrystals not valid";
0088 assert(0);
0089 }
0090 }
0091
0092 Basic2DGenericPFlowPositionCalc(const Basic2DGenericPFlowPositionCalc&) = delete;
0093 Basic2DGenericPFlowPositionCalc& operator=(const Basic2DGenericPFlowPositionCalc&) = delete;
0094
0095 void calculateAndSetPosition(reco::PFCluster&) override;
0096 void calculateAndSetPositions(reco::PFClusterCollection&) override;
0097
0098 private:
0099 const int _posCalcNCrystals;
0100 std::tuple<std::vector<int>, std::vector<int>, std::vector<float> > _logWeightDenom;
0101 const float _minAllowedNorm;
0102
0103 std::unique_ptr<CaloRecHitResolutionProvider> _timeResolutionCalcBarrel;
0104 std::unique_ptr<CaloRecHitResolutionProvider> _timeResolutionCalcEndcap;
0105
0106 void calculateAndSetPositionActual(reco::PFCluster&) const;
0107 };
0108
0109 DEFINE_EDM_PLUGIN(PFCPositionCalculatorFactory, Basic2DGenericPFlowPositionCalc, "Basic2DGenericPFlowPositionCalc");
0110
0111 namespace {
0112 inline bool isBarrel(int cell_layer) {
0113 return (cell_layer == PFLayer::HCAL_BARREL1 || cell_layer == PFLayer::HCAL_BARREL2 ||
0114 cell_layer == PFLayer::ECAL_BARREL);
0115 }
0116 }
0117
0118 void Basic2DGenericPFlowPositionCalc::calculateAndSetPosition(reco::PFCluster& cluster) {
0119 calculateAndSetPositionActual(cluster);
0120 }
0121
0122 void Basic2DGenericPFlowPositionCalc::calculateAndSetPositions(reco::PFClusterCollection& clusters) {
0123 for (reco::PFCluster& cluster : clusters) {
0124 calculateAndSetPositionActual(cluster);
0125 }
0126 }
0127
0128 void Basic2DGenericPFlowPositionCalc::calculateAndSetPositionActual(reco::PFCluster& cluster) const {
0129 if (!cluster.seed()) {
0130 throw cms::Exception("ClusterWithNoSeed") << " Found a cluster with no seed: " << cluster;
0131 }
0132 double cl_energy = 0;
0133 double cl_time = 0;
0134 double cl_timeweight = 0.0;
0135 double max_e = 0.0;
0136 PFLayer::Layer max_e_layer = PFLayer::NONE;
0137
0138
0139
0140
0141
0142
0143
0144 auto const recHitCollection =
0145 &(*cluster.recHitFractions()[0].recHitRef()) - cluster.recHitFractions()[0].recHitRef().key();
0146 auto nhits = cluster.recHitFractions().size();
0147 struct LHit {
0148 reco::PFRecHit const* hit;
0149 float energy;
0150 float fraction;
0151 };
0152 declareDynArray(LHit, nhits, hits);
0153 for (auto i = 0U; i < nhits; ++i) {
0154 auto const& hf = cluster.recHitFractions()[i];
0155 auto k = hf.recHitRef().key();
0156 auto p = recHitCollection + k;
0157 hits[i] = {p, (*p).energy(), float(hf.fraction())};
0158 }
0159
0160 bool resGiven = bool(_timeResolutionCalcBarrel) && bool(_timeResolutionCalcEndcap);
0161 LHit mySeed = {};
0162 for (auto const& rhf : hits) {
0163 const reco::PFRecHit& refhit = *rhf.hit;
0164 if (refhit.detId() == cluster.seed())
0165 mySeed = rhf;
0166 const auto rh_fraction = rhf.fraction;
0167 const auto rh_rawenergy = rhf.energy;
0168 const auto rh_energy = rh_rawenergy * rh_fraction;
0169 #ifdef PF_DEBUG
0170 if UNLIKELY (edm::isNotFinite(rh_energy)) {
0171 throw cms::Exception("PFClusterAlgo") << "rechit " << refhit.detId() << " has a NaN energy... "
0172 << "The input of the particle flow clustering seems to be corrupted.";
0173 }
0174 #endif
0175 cl_energy += rh_energy;
0176
0177 if (resGiven) {
0178 double res2 = 1.e-4;
0179 int cell_layer = (int)refhit.layer();
0180 res2 = isBarrel(cell_layer) ? 1. / _timeResolutionCalcBarrel->timeResolution2(rh_rawenergy)
0181 : 1. / _timeResolutionCalcEndcap->timeResolution2(rh_rawenergy);
0182 cl_time += rh_fraction * refhit.time() * res2;
0183 cl_timeweight += rh_fraction * res2;
0184 } else {
0185 const double rh_rawenergy2 = rh_rawenergy * rh_rawenergy;
0186 cl_timeweight += rh_rawenergy2 * rh_fraction;
0187 cl_time += rh_rawenergy2 * rh_fraction * refhit.time();
0188 }
0189
0190 if (rh_energy > max_e) {
0191 max_e = rh_energy;
0192 max_e_layer = refhit.layer();
0193 }
0194 }
0195
0196 cluster.setEnergy(cl_energy);
0197 cluster.setTime(cl_time / cl_timeweight);
0198 if (resGiven) {
0199 cluster.setTimeError(std::sqrt(1.0f / float(cl_timeweight)));
0200 }
0201 cluster.setLayer(max_e_layer);
0202
0203
0204 bool single_depth = true;
0205 int ref_depth = -1;
0206 double depth = 0.0;
0207 double position_norm = 0.0;
0208 double x(0.0), y(0.0), z(0.0);
0209 if (nullptr != mySeed.hit) {
0210 auto seedNeighbours = mySeed.hit->neighbours();
0211 switch (_posCalcNCrystals) {
0212 case 5:
0213 seedNeighbours = mySeed.hit->neighbours4();
0214 break;
0215 case 9:
0216 seedNeighbours = mySeed.hit->neighbours8();
0217 break;
0218 default:
0219 break;
0220 }
0221
0222 auto compute = [&](LHit const& rhf) {
0223 const reco::PFRecHit& refhit = *rhf.hit;
0224
0225 int cell_layer = (int)refhit.layer();
0226 float threshold = 0;
0227
0228 for (unsigned int j = 0; j < (std::get<2>(_logWeightDenom)).size(); ++j) {
0229
0230 int detectorEnum = std::get<0>(_logWeightDenom)[j];
0231 int depth = std::get<1>(_logWeightDenom)[j];
0232
0233 if ((cell_layer == PFLayer::HCAL_BARREL1 && detectorEnum == 1 && refhit.depth() == depth) ||
0234 (cell_layer == PFLayer::HCAL_ENDCAP && detectorEnum == 2 && refhit.depth() == depth) || detectorEnum == 0)
0235 threshold = std::get<2>(_logWeightDenom)[j];
0236 }
0237
0238 if (ref_depth < 0)
0239 ref_depth = refhit.depth();
0240 else if (refhit.depth() != ref_depth) {
0241
0242 single_depth = false;
0243 }
0244 const auto rh_energy = rhf.energy * rhf.fraction;
0245 const auto norm =
0246 (rhf.fraction < _minFractionInCalc ? 0.0f : std::max(0.0f, vdt::fast_logf(rh_energy * threshold)));
0247 const auto rhpos_xyz = refhit.position() * norm;
0248 x += rhpos_xyz.x();
0249 y += rhpos_xyz.y();
0250 z += rhpos_xyz.z();
0251 depth += refhit.depth() * norm;
0252 position_norm += norm;
0253 };
0254
0255 if (_posCalcNCrystals != -1)
0256 std::sort(hits.begin(), hits.end(), [](LHit const& a, LHit const& b) { return a.hit < b.hit; });
0257
0258 if (_posCalcNCrystals == -1)
0259 for (auto const& rhf : hits)
0260 compute(rhf);
0261 else {
0262 compute(mySeed);
0263
0264 unInitDynArray(reco::PFRecHit const*, seedNeighbours.size(), nei);
0265 for (auto k : seedNeighbours) {
0266 nei.push_back(recHitCollection + k);
0267 }
0268 std::sort(nei.begin(), nei.end());
0269 struct LHitLess {
0270 auto operator()(LHit const& a, reco::PFRecHit const* b) const { return a.hit < b; }
0271 auto operator()(reco::PFRecHit const* b, LHit const& a) const { return b < a.hit; }
0272 };
0273 std::set_intersection(
0274 hits.begin(), hits.end(), nei.begin(), nei.end(), boost::make_function_output_iterator(compute), LHitLess());
0275 }
0276 } else {
0277 throw cms::Exception("Basic2DGenerticPFlowPositionCalc")
0278 << "Cluster seed hit is null, something is wrong with PFlow RecHit!";
0279 }
0280
0281 if (position_norm < _minAllowedNorm) {
0282 edm::LogError("WeirdClusterNormalization") << "PFCluster too far from seeding cell: set position to (0,0,0).";
0283 cluster.setPosition(math::XYZPoint(0, 0, 0));
0284 cluster.calculatePositionREP();
0285 } else {
0286 const double norm_inverse = 1.0 / position_norm;
0287 x *= norm_inverse;
0288 y *= norm_inverse;
0289 z *= norm_inverse;
0290 if (single_depth)
0291 depth = ref_depth;
0292 else
0293 depth *= norm_inverse;
0294 cluster.setPosition(math::XYZPoint(x, y, z));
0295 cluster.setDepth(depth);
0296 cluster.calculatePositionREP();
0297 }
0298 }