CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
CTPPSTimingTrackRecognition.h
Go to the documentation of this file.
1 /****************************************************************************
2  *
3  * This is a part of CTPPS offline software.
4  * Authors:
5  * Laurent Forthomme (laurent.forthomme@cern.ch)
6  * Nicola Minafra (nicola.minafra@cern.ch)
7  * Mateusz Szpyrka (mateusz.szpyrka@cern.ch)
8  *
9  ****************************************************************************/
10 
11 #ifndef RecoPPS_Local_CTPPSTimingTrackRecognition
12 #define RecoPPS_Local_CTPPSTimingTrackRecognition
13 
16 
18 
20 
21 #include <vector>
22 #include <unordered_map>
23 #include <algorithm>
24 #include <cmath>
25 
30 template <typename TRACK_TYPE, typename HIT_TYPE>
32 public:
34  : threshold_(iConfig.getParameter<double>("threshold")),
35  thresholdFromMaximum_(iConfig.getParameter<double>("thresholdFromMaximum")),
36  resolution_(iConfig.getParameter<double>("resolution")),
37  sigma_(iConfig.getParameter<double>("sigma")),
38  tolerance_(iConfig.getParameter<double>("tolerance")),
39  pixelEfficiencyFunction_(iConfig.getParameter<std::string>("pixelEfficiencyFunction")) {
41  throw cms::Exception("CTPPSTimingTrackRecognition")
42  << "Invalid number of parameters to the pixel efficiency function! "
44  }
45  virtual ~CTPPSTimingTrackRecognition() = default;
46 
47  //--- class API
48 
50  inline virtual void clear() { hitVectorMap_.clear(); }
52  virtual void addHit(const HIT_TYPE& recHit) = 0;
55 
56 protected:
57  // Algorithm parameters:
58  const float threshold_;
59  const float thresholdFromMaximum_;
60  const float resolution_;
61  const float sigma_;
62  const float tolerance_;
64 
65  typedef std::vector<TRACK_TYPE> TrackVector;
66  typedef std::vector<HIT_TYPE> HitVector;
67  typedef std::unordered_map<int, HitVector> HitVectorMap;
68 
71 
76  };
78  struct SpatialRange {
79  float xBegin, xEnd;
80  float yBegin, yEnd;
81  float zBegin, zEnd;
82  };
83 
94  const DimensionParameters& param,
95  float (*getHitCenter)(const HIT_TYPE&),
96  float (*getHitRangeWidth)(const HIT_TYPE&),
97  void (*setTrackCenter)(TRACK_TYPE&, float),
98  void (*setTrackSigma)(TRACK_TYPE&, float),
100 
111  bool timeEval(const HitVector& hits, float& meanTime, float& timeSigma) const;
112 };
113 
114 /****************************************************************************
115  * Implementation
116  ****************************************************************************/
117 
118 template <class TRACK_TYPE, class HIT_TYPE>
120  const HitVector& hits,
121  const DimensionParameters& param,
122  float (*getHitCenter)(const HIT_TYPE&),
123  float (*getHitRangeWidth)(const HIT_TYPE&),
124  void (*setTrackCenter)(TRACK_TYPE&, float),
125  void (*setTrackSigma)(TRACK_TYPE&, float),
126  TrackVector& result) {
127  int numberOfTracks = 0;
128  const float invResolution = 1. / resolution_;
129  const float profileRangeMargin = sigma_ * 3.;
130  const float profileRangeBegin = param.rangeBegin - profileRangeMargin;
131  const float profileRangeEnd = param.rangeEnd + profileRangeMargin;
132 
133  std::vector<float> hitProfile((profileRangeEnd - profileRangeBegin) * invResolution + 1, 0.);
134  // extra component to make sure that the profile drops below the threshold at range's end
135  *hitProfile.rbegin() = -1.f;
136 
137  // Creates hit profile
138  for (auto const& hit : hits) {
139  const float center = getHitCenter(hit), rangeWidth = getHitRangeWidth(hit);
140  std::vector<double> params{center, rangeWidth, sigma_};
141  for (unsigned int i = 0; i < hitProfile.size(); ++i)
142  hitProfile[i] +=
143  pixelEfficiencyFunction_.evaluate(std::vector<double>{profileRangeBegin + i * resolution_}, params);
144  }
145 
146  bool underThreshold = true;
147  float rangeMaximum = -1.0f;
148  bool trackRangeFound = false;
149  int trackRangeBegin = 0, trackRangeEnd = 0;
150 
151  // Searches for tracks in the hit profile
152  for (unsigned int i = 0; i < hitProfile.size(); i++) {
153  if (hitProfile[i] > rangeMaximum)
154  rangeMaximum = hitProfile[i];
155 
156  // Going above the threshold
157  if (underThreshold && hitProfile[i] > threshold_) {
158  underThreshold = false;
159  trackRangeBegin = i;
160  rangeMaximum = hitProfile[i];
161  }
162 
163  // Going under the threshold
164  else if (!underThreshold && hitProfile[i] <= threshold_) {
165  underThreshold = true;
166  trackRangeEnd = i;
167  trackRangeFound = true;
168  }
169 
170  // Finds all tracks within the track range
171  if (trackRangeFound) {
172  float trackThreshold = rangeMaximum - thresholdFromMaximum_;
173  int trackBegin;
174  bool underTrackThreshold = true;
175 
176  for (int j = trackRangeBegin; j <= trackRangeEnd; j++) {
177  if (underTrackThreshold && hitProfile[j] > trackThreshold) {
178  underTrackThreshold = false;
179  trackBegin = j;
180  } else if (!underTrackThreshold && hitProfile[j] <= trackThreshold) {
181  underTrackThreshold = true;
182  TRACK_TYPE track;
183  float leftMargin = profileRangeBegin + resolution_ * trackBegin;
184  float rightMargin = profileRangeBegin + resolution_ * j;
185  setTrackCenter(track, 0.5f * (leftMargin + rightMargin));
186  setTrackSigma(track, 0.5f * (rightMargin - leftMargin));
187  result.push_back(track);
188  numberOfTracks++;
189  }
190  }
191  trackRangeFound = false;
192  }
193  }
194 }
195 
196 template <class TRACK_TYPE, class HIT_TYPE>
199  bool initialized = false;
200  SpatialRange result = {};
201 
202  for (const auto& hit : hits) {
203  const float xBegin = hit.x() - 0.5f * hit.xWidth(), xEnd = hit.x() + 0.5f * hit.xWidth();
204  const float yBegin = hit.y() - 0.5f * hit.yWidth(), yEnd = hit.y() + 0.5f * hit.yWidth();
205  const float zBegin = hit.z() - 0.5f * hit.zWidth(), zEnd = hit.z() + 0.5f * hit.zWidth();
206 
207  if (!initialized) {
208  result.xBegin = xBegin;
209  result.xEnd = xEnd;
210  result.yBegin = yBegin;
211  result.yEnd = yEnd;
212  result.zBegin = zBegin;
213  result.zEnd = zEnd;
214  initialized = true;
215  continue;
216  }
217  result.xBegin = std::min(result.xBegin, xBegin);
218  result.xEnd = std::max(result.xEnd, xEnd);
219  result.yBegin = std::min(result.yBegin, yBegin);
220  result.yEnd = std::max(result.yEnd, yEnd);
221  result.zBegin = std::min(result.zBegin, zBegin);
222  result.zEnd = std::max(result.zEnd, zEnd);
223  }
224 
225  return result;
226 }
227 
228 template <class TRACK_TYPE, class HIT_TYPE>
230  float& mean_time,
231  float& time_sigma) const {
232  float mean_num = 0.f, mean_denom = 0.f;
233  bool valid_hits = false;
234  for (const auto& hit : hits) {
235  if (hit.tPrecision() <= 0.)
236  continue;
237  const float weight = std::pow(hit.tPrecision(), -2);
238  mean_num += weight * hit.time();
239  mean_denom += weight;
240  valid_hits = true; // at least one valid hit to account for
241  }
242  mean_time = valid_hits ? (mean_num / mean_denom) : 0.f;
243  time_sigma = valid_hits ? std::sqrt(1.f / mean_denom) : 0.f;
244  return valid_hits;
245 }
246 
247 #endif
SpatialRange getHitSpatialRange(const HitVector &hits)
double evaluate(V const &iVariables, P const &iParameters) const
Structure representing a 3D range in space.
void producePartialTracks(const HitVector &hits, const DimensionParameters &param, float(*getHitCenter)(const HIT_TYPE &), float(*getHitRangeWidth)(const HIT_TYPE &), void(*setTrackCenter)(TRACK_TYPE &, float), void(*setTrackSigma)(TRACK_TYPE &, float), TrackVector &result)
virtual ~CTPPSTimingTrackRecognition()=default
Definition: weight.py:1
virtual void clear()
Reset internal state of a class instance.
CTPPSTimingTrackRecognition(const edm::ParameterSet &iConfig)
std::unordered_map< int, HitVector > HitVectorMap
T sqrt(T t)
Definition: SSEVec.h:19
double f[11][100]
unsigned int numberOfParameters() const
reco::FormulaEvaluator pixelEfficiencyFunction_
std::vector< TRACK_TYPE > TrackVector
HitVectorMap hitVectorMap_
RecHit vectors that should be processed separately while reconstructing tracks.
virtual void addHit(const HIT_TYPE &recHit)=0
Add new hit to the set from which the tracks are reconstructed.
bool timeEval(const HitVector &hits, float &meanTime, float &timeSigma) const
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
virtual int produceTracks(edm::DetSet< TRACK_TYPE > &tracks)=0
Produce a collection of tracks, given its hits collection.