File indexing completed on 2023-10-25 10:01:54
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011 #ifndef RecoPPS_Local_CTPPSTimingTrackRecognition
0012 #define RecoPPS_Local_CTPPSTimingTrackRecognition
0013
0014 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0015 #include "FWCore/Utilities/interface/Exception.h"
0016
0017 #include "CommonTools/Utils/interface/FormulaEvaluator.h"
0018
0019 #include "DataFormats/Common/interface/DetSet.h"
0020
0021 #include <vector>
0022 #include <unordered_map>
0023 #include <algorithm>
0024 #include <cmath>
0025
0026
0027
0028
0029
0030 template <typename TRACK_TYPE, typename HIT_TYPE>
0031 class CTPPSTimingTrackRecognition {
0032 public:
0033 inline CTPPSTimingTrackRecognition(const edm::ParameterSet& iConfig)
0034 : threshold_(iConfig.getParameter<double>("threshold")),
0035 thresholdFromMaximum_(iConfig.getParameter<double>("thresholdFromMaximum")),
0036 resolution_(iConfig.getParameter<double>("resolution")),
0037 sigma_(iConfig.getParameter<double>("sigma")),
0038 tolerance_(iConfig.getParameter<double>("tolerance")),
0039 pixelEfficiencyFunction_(iConfig.getParameter<std::string>("pixelEfficiencyFunction")) {
0040 if (pixelEfficiencyFunction_.numberOfParameters() != 3)
0041 throw cms::Exception("CTPPSTimingTrackRecognition")
0042 << "Invalid number of parameters to the pixel efficiency function! "
0043 << pixelEfficiencyFunction_.numberOfParameters() << " != 3.";
0044 }
0045 virtual ~CTPPSTimingTrackRecognition() = default;
0046
0047
0048
0049
0050 inline virtual void clear() { hitVectorMap_.clear(); }
0051
0052 virtual void addHit(const HIT_TYPE& recHit) = 0;
0053
0054 virtual int produceTracks(edm::DetSet<TRACK_TYPE>& tracks) = 0;
0055
0056 protected:
0057
0058 const float threshold_;
0059 const float thresholdFromMaximum_;
0060 const float resolution_;
0061 const float sigma_;
0062 const float tolerance_;
0063 reco::FormulaEvaluator pixelEfficiencyFunction_;
0064
0065 typedef std::vector<TRACK_TYPE> TrackVector;
0066 typedef std::vector<HIT_TYPE> HitVector;
0067 typedef std::unordered_map<int, HitVector> HitVectorMap;
0068
0069
0070 HitVectorMap hitVectorMap_;
0071
0072
0073
0074 struct DimensionParameters {
0075 float rangeBegin, rangeEnd;
0076 };
0077
0078 struct SpatialRange {
0079 float xBegin, xEnd;
0080 float yBegin, yEnd;
0081 float zBegin, zEnd;
0082 };
0083
0084
0085
0086
0087
0088
0089
0090
0091
0092
0093 void producePartialTracks(const HitVector& hits,
0094 const DimensionParameters& param,
0095 float (*getHitCenter)(const HIT_TYPE&),
0096 float (*getHitRangeWidth)(const HIT_TYPE&),
0097 void (*setTrackCenter)(TRACK_TYPE&, float),
0098 void (*setTrackSigma)(TRACK_TYPE&, float),
0099 TrackVector& result);
0100
0101
0102
0103
0104 SpatialRange getHitSpatialRange(const HitVector& hits);
0105
0106
0107
0108
0109
0110
0111 bool timeEval(const HitVector& hits, float& meanTime, float& timeSigma) const;
0112 };
0113
0114
0115
0116
0117
0118 template <class TRACK_TYPE, class HIT_TYPE>
0119 inline void CTPPSTimingTrackRecognition<TRACK_TYPE, HIT_TYPE>::producePartialTracks(
0120 const HitVector& hits,
0121 const DimensionParameters& param,
0122 float (*getHitCenter)(const HIT_TYPE&),
0123 float (*getHitRangeWidth)(const HIT_TYPE&),
0124 void (*setTrackCenter)(TRACK_TYPE&, float),
0125 void (*setTrackSigma)(TRACK_TYPE&, float),
0126 TrackVector& result) {
0127 const float invResolution = 1. / resolution_;
0128 const float profileRangeMargin = sigma_ * 3.;
0129 const float profileRangeBegin = param.rangeBegin - profileRangeMargin;
0130 const float profileRangeEnd = param.rangeEnd + profileRangeMargin;
0131
0132 std::vector<float> hitProfile((profileRangeEnd - profileRangeBegin) * invResolution + 1, 0.);
0133
0134 *hitProfile.rbegin() = -1.f;
0135
0136
0137 for (auto const& hit : hits) {
0138 const float center = getHitCenter(hit), rangeWidth = getHitRangeWidth(hit);
0139 std::vector<double> params{center, rangeWidth, sigma_};
0140 for (unsigned int i = 0; i < hitProfile.size(); ++i)
0141 hitProfile[i] +=
0142 pixelEfficiencyFunction_.evaluate(std::vector<double>{profileRangeBegin + i * resolution_}, params);
0143 }
0144
0145 bool underThreshold = true;
0146 float rangeMaximum = -1.0f;
0147 bool trackRangeFound = false;
0148 int trackRangeBegin = 0, trackRangeEnd = 0;
0149
0150
0151 for (unsigned int i = 0; i < hitProfile.size(); i++) {
0152 if (hitProfile[i] > rangeMaximum)
0153 rangeMaximum = hitProfile[i];
0154
0155
0156 if (underThreshold && hitProfile[i] > threshold_) {
0157 underThreshold = false;
0158 trackRangeBegin = i;
0159 rangeMaximum = hitProfile[i];
0160 }
0161
0162
0163 else if (!underThreshold && hitProfile[i] <= threshold_) {
0164 underThreshold = true;
0165 trackRangeEnd = i;
0166 trackRangeFound = true;
0167 }
0168
0169
0170 if (trackRangeFound) {
0171 float trackThreshold = rangeMaximum - thresholdFromMaximum_;
0172 int trackBegin;
0173 bool underTrackThreshold = true;
0174
0175 for (int j = trackRangeBegin; j <= trackRangeEnd; j++) {
0176 if (underTrackThreshold && hitProfile[j] > trackThreshold) {
0177 underTrackThreshold = false;
0178 trackBegin = j;
0179 } else if (!underTrackThreshold && hitProfile[j] <= trackThreshold) {
0180 underTrackThreshold = true;
0181 TRACK_TYPE track;
0182 float leftMargin = profileRangeBegin + resolution_ * trackBegin;
0183 float rightMargin = profileRangeBegin + resolution_ * j;
0184 setTrackCenter(track, 0.5f * (leftMargin + rightMargin));
0185 setTrackSigma(track, 0.5f * (rightMargin - leftMargin));
0186 result.push_back(track);
0187 }
0188 }
0189 trackRangeFound = false;
0190 }
0191 }
0192 }
0193
0194 template <class TRACK_TYPE, class HIT_TYPE>
0195 inline typename CTPPSTimingTrackRecognition<TRACK_TYPE, HIT_TYPE>::SpatialRange
0196 CTPPSTimingTrackRecognition<TRACK_TYPE, HIT_TYPE>::getHitSpatialRange(const HitVector& hits) {
0197 bool initialized = false;
0198 SpatialRange result = {};
0199
0200 for (const auto& hit : hits) {
0201 const float xBegin = hit.x() - 0.5f * hit.xWidth(), xEnd = hit.x() + 0.5f * hit.xWidth();
0202 const float yBegin = hit.y() - 0.5f * hit.yWidth(), yEnd = hit.y() + 0.5f * hit.yWidth();
0203 const float zBegin = hit.z() - 0.5f * hit.zWidth(), zEnd = hit.z() + 0.5f * hit.zWidth();
0204
0205 if (!initialized) {
0206 result.xBegin = xBegin;
0207 result.xEnd = xEnd;
0208 result.yBegin = yBegin;
0209 result.yEnd = yEnd;
0210 result.zBegin = zBegin;
0211 result.zEnd = zEnd;
0212 initialized = true;
0213 continue;
0214 }
0215 result.xBegin = std::min(result.xBegin, xBegin);
0216 result.xEnd = std::max(result.xEnd, xEnd);
0217 result.yBegin = std::min(result.yBegin, yBegin);
0218 result.yEnd = std::max(result.yEnd, yEnd);
0219 result.zBegin = std::min(result.zBegin, zBegin);
0220 result.zEnd = std::max(result.zEnd, zEnd);
0221 }
0222
0223 return result;
0224 }
0225
0226 template <class TRACK_TYPE, class HIT_TYPE>
0227 inline bool CTPPSTimingTrackRecognition<TRACK_TYPE, HIT_TYPE>::timeEval(const HitVector& hits,
0228 float& mean_time,
0229 float& time_sigma) const {
0230 float mean_num = 0.f, mean_denom = 0.f;
0231 bool valid_hits = false;
0232 for (const auto& hit : hits) {
0233 if (hit.tPrecision() <= 0.)
0234 continue;
0235 const float weight = std::pow(hit.tPrecision(), -2);
0236 mean_num += weight * hit.time();
0237 mean_denom += weight;
0238 valid_hits = true;
0239 }
0240 mean_time = valid_hits ? (mean_num / mean_denom) : 0.f;
0241 time_sigma = valid_hits ? std::sqrt(1.f / mean_denom) : 0.f;
0242 return valid_hits;
0243 }
0244
0245 #endif