PFClusterProducer/plugins/Basic2DGenericPFlowPositionCalc.cc

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