plugins/SimMappers/RealisticHitToClusterAssociator.h

0001 #ifndef __RecoParticleFlow_PFClusterProducer_RealisticHitToClusterAssociator_H__
0002 #define __RecoParticleFlow_PFClusterProducer_RealisticHitToClusterAssociator_H__
0003 /////////////////////////
0004 // Author: Felice Pantaleo
0005 // Date:   30/06/2017
0006 // Email: felice@cern.ch
0007 /////////////////////////
0008 #include <vector>
0009 #include <cmath>
0010 #include <unordered_map>
0011 #include "RealisticCluster.h"
0012
0013 class RealisticHitToClusterAssociator {
0014   using Hit3DPosition = std::array<float, 3>;
0015
0016 public:
0017   struct RealisticHit {
0018     struct HitToCluster {
0019       unsigned int simClusterId_;
0020       float mcEnergyFraction_;
0021       float distanceFromMaxHit_;
0022       float realisticEnergyFraction_;
0023     };
0024
0025     Hit3DPosition hitPosition_;
0026     float totalEnergy_;
0027     unsigned int layerId_;
0028     std::vector<HitToCluster> hitToCluster_;
0029   };
0030
0031   void init(std::size_t numberOfHits, std::size_t numberOfSimClusters, std::size_t numberOfLayers) {
0032     realisticHits_.resize(numberOfHits);
0033     realisticSimClusters_.resize(numberOfSimClusters);
0034     for (auto& sc : realisticSimClusters_)
0035       sc.setLayersNum(numberOfLayers);
0036   }
0037
0038   void insertHitPosition(float x, float y, float z, unsigned int hitIndex) {
0039     realisticHits_[hitIndex].hitPosition_ = {{x, y, z}};
0040   }
0041
0042   void insertLayerId(unsigned int layerId, unsigned int hitIndex) { realisticHits_[hitIndex].layerId_ = layerId; }
0043
0044   void insertHitEnergy(float energy, unsigned int hitIndex) { realisticHits_[hitIndex].totalEnergy_ = energy; }
0045
0046   void insertSimClusterIdAndFraction(unsigned int scIdx, float fraction, unsigned int hitIndex, float associatedEnergy) {
0047     realisticHits_[hitIndex].hitToCluster_.emplace_back(RealisticHit::HitToCluster{scIdx, fraction, 0.f, 0.f});
0048     realisticSimClusters_[scIdx].setMaxEnergyHit(
0049         realisticHits_[hitIndex].layerId_, associatedEnergy, realisticHits_[hitIndex].hitPosition_);
0050   }
0051
0052   float XYdistanceFromMaxHit(unsigned int hitId, unsigned int clusterId) {
0053     auto l = realisticHits_[hitId].layerId_;
0054     const auto& maxHitPosition = realisticSimClusters_[clusterId].getMaxEnergyPosition(l);
0055     float distanceSquared = std::pow((realisticHits_[hitId].hitPosition_[0] - maxHitPosition[0]), 2) +
0056                             std::pow((realisticHits_[hitId].hitPosition_[1] - maxHitPosition[1]), 2);
0057     return std::sqrt(distanceSquared);
0058   }
0059
0060   float XYdistanceFromPointOnSameLayer(unsigned int hitId, const Hit3DPosition& point) {
0061     float distanceSquared = std::pow((realisticHits_[hitId].hitPosition_[0] - point[0]), 2) +
0062                             std::pow((realisticHits_[hitId].hitPosition_[1] - point[1]), 2);
0063     return std::sqrt(distanceSquared);
0064   }
0065
0066   void computeAssociation(float exclusiveFraction,
0067                           bool useMCFractionsForExclEnergy,
0068                           unsigned int fhOffset,
0069                           unsigned int bhOffset) {
0070     //if more than exclusiveFraction of a hit's energy belongs to a cluster, that rechit is not counted as shared
0071     unsigned int numberOfHits = realisticHits_.size();
0072     std::vector<float> partialEnergies;
0073     for (unsigned int hitId = 0; hitId < numberOfHits; ++hitId) {
0074       partialEnergies.clear();
0075       std::vector<unsigned int> removeAssociation;
0076       auto& realisticHit = realisticHits_[hitId];
0077       unsigned int numberOfClusters = realisticHit.hitToCluster_.size();
0078       if (numberOfClusters == 1) {
0079         unsigned int simClusterId = realisticHit.hitToCluster_[0].simClusterId_;
0080         float assignedFraction = 1.f;
0081         realisticHit.hitToCluster_[0].realisticEnergyFraction_ = assignedFraction;
0082         float assignedEnergy = realisticHit.totalEnergy_;
0083         realisticSimClusters_[simClusterId].increaseEnergy(assignedEnergy);
0084         realisticSimClusters_[simClusterId].addHitAndFraction(hitId, assignedFraction);
0085         realisticSimClusters_[simClusterId].increaseExclusiveEnergy(assignedEnergy);
0086       } else {
0087         partialEnergies.resize(numberOfClusters, 0.f);
0088         unsigned int layer = realisticHit.layerId_;
0089         float sumE = 0.f;
0090         float energyDecayLength = getDecayLength(layer, fhOffset, bhOffset);
0091         for (unsigned int clId = 0; clId < numberOfClusters; ++clId) {
0092           auto simClusterId = realisticHit.hitToCluster_[clId].simClusterId_;
0093           realisticHit.hitToCluster_[clId].distanceFromMaxHit_ = XYdistanceFromMaxHit(hitId, simClusterId);
0094           // partial energy is computed based on the distance from the maximum energy hit and its energy
0095           // partial energy is only needed to compute a fraction and it's not the energy assigned to the cluster
0096           auto maxEnergyOnLayer = realisticSimClusters_[simClusterId].getMaxEnergy(layer);
0097           if (maxEnergyOnLayer > 0.f) {
0098             partialEnergies[clId] =
0099                 maxEnergyOnLayer * std::exp(-realisticHit.hitToCluster_[clId].distanceFromMaxHit_ / energyDecayLength);
0100           }
0101           sumE += partialEnergies[clId];
0102         }
0103         if (sumE > 0.f) {
0104           float invSumE = 1.f / sumE;
0105           for (unsigned int clId = 0; clId < numberOfClusters; ++clId) {
0106             unsigned int simClusterIndex = realisticHit.hitToCluster_[clId].simClusterId_;
0107             float assignedFraction = partialEnergies[clId] * invSumE;
0108             if (assignedFraction > 1e-3) {
0109               realisticHit.hitToCluster_[clId].realisticEnergyFraction_ = assignedFraction;
0110               float assignedEnergy = assignedFraction * realisticHit.totalEnergy_;
0111               realisticSimClusters_[simClusterIndex].increaseEnergy(assignedEnergy);
0112               realisticSimClusters_[simClusterIndex].addHitAndFraction(hitId, assignedFraction);
0113               // if the hits energy belongs for more than exclusiveFraction to a cluster, also the cluster's
0114               // exclusive energy is increased. The exclusive energy will be needed to evaluate if
0115               // a realistic cluster will be invisible, i.e. absorbed by other clusters
0116               if ((useMCFractionsForExclEnergy and
0117                    realisticHit.hitToCluster_[clId].mcEnergyFraction_ > exclusiveFraction) or
0118                   (!useMCFractionsForExclEnergy and assignedFraction > exclusiveFraction)) {
0119                 realisticSimClusters_[simClusterIndex].increaseExclusiveEnergy(assignedEnergy);
0120               }
0121             } else {
0122               removeAssociation.push_back(simClusterIndex);
0123             }
0124           }
0125         }
0126       }
0127
0128       while (!removeAssociation.empty()) {
0129         auto clusterToRemove = removeAssociation.back();
0130         removeAssociation.pop_back();
0131
0132         realisticHit.hitToCluster_.erase(std::remove_if(realisticHit.hitToCluster_.begin(),
0133                                                         realisticHit.hitToCluster_.end(),
0134                                                         [clusterToRemove](const RealisticHit::HitToCluster& x) {
0135                                                           return x.simClusterId_ == clusterToRemove;
0136                                                         }),
0137                                          realisticHit.hitToCluster_.end());
0138       }
0139     }
0140   }
0141
0142   void findAndMergeInvisibleClusters(float invisibleFraction, float exclusiveFraction) {
0143     unsigned int numberOfRealSimClusters = realisticSimClusters_.size();
0144
0145     for (unsigned int clId = 0; clId < numberOfRealSimClusters; ++clId) {
0146       if (realisticSimClusters_[clId].getExclusiveEnergyFraction() < invisibleFraction) {
0147         realisticSimClusters_[clId].setVisible(false);
0148         auto& hAndF = realisticSimClusters_[clId].hitsIdsAndFractions();
0149         std::unordered_map<unsigned int, float> energyInNeighbors;
0150         float totalSharedEnergy = 0.f;
0151
0152         for (auto& elt : hAndF) {
0153           unsigned int hitId = elt.first;
0154           float fraction = elt.second;
0155           auto& realisticHit = realisticHits_[hitId];
0156
0157           if (realisticHit.hitToCluster_.size() > 1 && fraction < 1.f) {
0158             float correction = 1.f - fraction;
0159             unsigned int numberOfClusters = realisticHit.hitToCluster_.size();
0160             int clusterToRemove = -1;
0161             for (unsigned int i = 0; i < numberOfClusters; ++i) {
0162               auto simClusterIndex = realisticHit.hitToCluster_[i].simClusterId_;
0163               if (simClusterIndex == clId) {
0164                 clusterToRemove = i;
0165               } else if (realisticSimClusters_[simClusterIndex].isVisible()) {
0166                 float oldFraction = realisticHit.hitToCluster_[i].realisticEnergyFraction_;
0167                 float newFraction = oldFraction / correction;
0168                 float oldEnergy = oldFraction * realisticHit.totalEnergy_;
0169                 float newEnergy = newFraction * realisticHit.totalEnergy_;
0170                 float sharedEnergy = newEnergy - oldEnergy;
0171                 energyInNeighbors[simClusterIndex] += sharedEnergy;
0172                 totalSharedEnergy += sharedEnergy;
0173                 realisticSimClusters_[simClusterIndex].increaseEnergy(sharedEnergy);
0174                 realisticSimClusters_[simClusterIndex].modifyFractionForHitId(newFraction, hitId);
0175                 realisticHit.hitToCluster_[i].realisticEnergyFraction_ = newFraction;
0176                 if (newFraction > exclusiveFraction) {
0177                   realisticSimClusters_[simClusterIndex].increaseExclusiveEnergy(sharedEnergy);
0178                   if (oldFraction <= exclusiveFraction) {
0179                     realisticSimClusters_[simClusterIndex].increaseExclusiveEnergy(oldEnergy);
0180                   }
0181                 }
0182               }
0183             }
0184             realisticSimClusters_[realisticHit.hitToCluster_[clusterToRemove].simClusterId_].modifyFractionForHitId(
0185                 0.f, hitId);
0186             realisticHit.hitToCluster_.erase(realisticHit.hitToCluster_.begin() + clusterToRemove);
0187           }
0188         }
0189
0190         for (auto& elt : hAndF) {
0191           unsigned int hitId = elt.first;
0192           auto& realisticHit = realisticHits_[hitId];
0193           if (realisticHit.hitToCluster_.size() == 1 and realisticHit.hitToCluster_[0].simClusterId_ == clId and
0194               totalSharedEnergy > 0.f) {
0195             for (auto& pair : energyInNeighbors) {
0196               // hits that belonged completely to the absorbed cluster are redistributed
0197               // based on the fraction of energy shared in the shared hits
0198               float sharedFraction = pair.second / totalSharedEnergy;
0199               if (sharedFraction > 1e-6) {
0200                 float assignedEnergy = realisticHit.totalEnergy_ * sharedFraction;
0201                 realisticSimClusters_[pair.first].increaseEnergy(assignedEnergy);
0202                 realisticSimClusters_[pair.first].addHitAndFraction(hitId, sharedFraction);
0203                 realisticHit.hitToCluster_.emplace_back(
0204                     RealisticHit::HitToCluster{pair.first, 0.f, -1.f, sharedFraction});
0205                 if (sharedFraction > exclusiveFraction)
0206                   realisticSimClusters_[pair.first].increaseExclusiveEnergy(assignedEnergy);
0207               }
0208             }
0209           }
0210         }
0211       }
0212     }
0213   }
0214
0215   void findCentersOfGravity() {
0216     for (auto& cluster : realisticSimClusters_) {
0217       if (cluster.isVisible()) {
0218         unsigned int layersNum = cluster.getLayersNum();
0219         std::vector<float> totalEnergyPerLayer(layersNum, 0.f);
0220         std::vector<float> xEnergyPerLayer(layersNum, 0.f);
0221         std::vector<float> yEnergyPerLayer(layersNum, 0.f);
0222         std::vector<float> zPositionPerLayer(layersNum, 0.f);
0223         const auto& hAndF = cluster.hitsIdsAndFractions();
0224         for (auto& elt : hAndF) {
0225           auto hitId = elt.first;
0226           auto fraction = elt.second;
0227           const auto& hit = realisticHits_[hitId];
0228           const auto& hitPos = hit.hitPosition_;
0229           auto layerId = hit.layerId_;
0230           auto hitEinCluster = hit.totalEnergy_ * fraction;
0231           totalEnergyPerLayer[layerId] += hitEinCluster;
0232           xEnergyPerLayer[layerId] += hitPos[0] * hitEinCluster;
0233           yEnergyPerLayer[layerId] += hitPos[1] * hitEinCluster;
0234           zPositionPerLayer[layerId] = hitPos[2];
0235         }
0236         Hit3DPosition centerOfGravity;
0237         for (unsigned int layerId = 0; layerId < layersNum; layerId++) {
0238           auto energyOnLayer = totalEnergyPerLayer[layerId];
0239           if (energyOnLayer > 0.f) {
0240             centerOfGravity = {{xEnergyPerLayer[layerId] / energyOnLayer,
0241                                 yEnergyPerLayer[layerId] / energyOnLayer,
0242                                 zPositionPerLayer[layerId]}};
0243             cluster.setCenterOfGravity(layerId, centerOfGravity);
0244           }
0245         }
0246       }
0247     }
0248   }
0249
0250   void filterHitsByDistance(float maxDistance) {
0251     for (auto& cluster : realisticSimClusters_) {
0252       if (cluster.isVisible()) {
0253         auto& hAndF = cluster.hitsIdsAndFractions();
0254         for (unsigned int i = 0; i < hAndF.size(); ++i) {
0255           auto hitId = hAndF[i].first;
0256           const auto& hit = realisticHits_[hitId];
0257           auto layerId = hit.layerId_;
0258           if (XYdistanceFromPointOnSameLayer(hitId, cluster.getCenterOfGravity(layerId)) > maxDistance) {
0259             cluster.increaseEnergy(-hit.totalEnergy_ * hAndF[i].second);
0260             cluster.modifyFractionByIndex(0.f, i);
0261           }
0262         }
0263       }
0264     }
0265   }
0266
0267   const std::vector<RealisticCluster>& realisticClusters() const { return realisticSimClusters_; }
0268
0269 private:
0270   static float getDecayLength(unsigned int layer, unsigned int fhOffset, unsigned int bhOffset) {
0271     constexpr float eeDecayLengthInLayer = 2.f;
0272     constexpr float fhDecayLengthInLayer = 1.5f;
0273     constexpr float bhDecayLengthInLayer = 1.f;
0274
0275     if (layer <= fhOffset)
0276       return eeDecayLengthInLayer;
0277     else if (layer > fhOffset && layer <= bhOffset)
0278       return fhDecayLengthInLayer;
0279     else
0280       return bhDecayLengthInLayer;
0281   }
0282
0283   // the vector of the Realistic SimClusters
0284   std::vector<RealisticCluster> realisticSimClusters_;
0285   // the vector of the Realistic SimClusters
0286   std::vector<RealisticHit> realisticHits_;
0287 };
0288
0289 #endif