Local/interface/CTPPSTimingTrackRecognition.h

0001 /****************************************************************************
0002  *
0003  * This is a part of CTPPS offline software.
0004  * Authors:
0005  *   Laurent Forthomme (laurent.forthomme@cern.ch)
0006  *   Nicola Minafra (nicola.minafra@cern.ch)
0007  *   Mateusz Szpyrka (mateusz.szpyrka@cern.ch)
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  * Class intended to perform general CTPPS timing detectors track recognition,
0028  * as well as construction of specialized classes (for now CTPPSDiamond and TotemTiming local tracks).
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   //--- class API
0048
0049   /// Reset internal state of a class instance.
0050   inline virtual void clear() { hitVectorMap_.clear(); }
0051   /// Add new hit to the set from which the tracks are reconstructed.
0052   virtual void addHit(const HIT_TYPE& recHit) = 0;
0053   /// Produce a collection of tracks, given its hits collection
0054   virtual int produceTracks(edm::DetSet<TRACK_TYPE>& tracks) = 0;
0055
0056 protected:
0057   // Algorithm parameters:
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   /// RecHit vectors that should be processed separately while reconstructing tracks.
0070   HitVectorMap hitVectorMap_;
0071
0072   /// Structure representing parameters set for single dimension.
0073   /// Intended to use when producing partial tracks.
0074   struct DimensionParameters {
0075     float rangeBegin, rangeEnd;
0076   };
0077   /// Structure representing a 3D range in space.
0078   struct SpatialRange {
0079     float xBegin, xEnd;
0080     float yBegin, yEnd;
0081     float zBegin, zEnd;
0082   };
0083
0084   /** Produce all partial tracks from given set with regard to single dimension.
0085      * \param[in] hits vector of hits from which the tracks are created
0086      * \param[in] param describe all parameters used by 1D track recognition algorithm
0087      * \param[in] getHitCenter function extracting hit's center in the dimension that the partial tracks relate to
0088      * \param[in] getHitRangeWidth analogue to \a getHitCenter, but extracts hit's width in specific dimension
0089      * \param[in] setTrackCenter function used to set track's position in considered dimension
0090      * \param[in] setTrackSigma function used to set track's sigma in considered dimension
0091      * \param[out] result vector to which produced tracks are appended
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   /** Retrieve the bounds of a 3D range in which all hits from given collection are contained.
0102      * \param[in] hits hits collection to retrieve the range from
0103      */
0104   SpatialRange getHitSpatialRange(const HitVector& hits);
0105   /** Evaluate the time + associated uncertainty for a given track
0106      * \note General remarks:
0107      * - track's time = weighted mean of all hit times with time precision as weight,
0108      * - track's time sigma = uncertainty of the weighted mean
0109      * - hit is ignored if the time precision is equal to 0
0110      */
0111   bool timeEval(const HitVector& hits, float& meanTime, float& timeSigma) const;
0112 };
0113
0114 /****************************************************************************
0115  * Implementation
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   // extra component to make sure that the profile drops below the threshold at range's end
0134   *hitProfile.rbegin() = -1.f;
0135
0136   // Creates hit profile
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   // Searches for tracks in the hit profile
0151   for (unsigned int i = 0; i < hitProfile.size(); i++) {
0152     if (hitProfile[i] > rangeMaximum)
0153       rangeMaximum = hitProfile[i];
0154
0155     // Going above the threshold
0156     if (underThreshold && hitProfile[i] > threshold_) {
0157       underThreshold = false;
0158       trackRangeBegin = i;
0159       rangeMaximum = hitProfile[i];
0160     }
0161
0162     // Going under the threshold
0163     else if (!underThreshold && hitProfile[i] <= threshold_) {
0164       underThreshold = true;
0165       trackRangeEnd = i;
0166       trackRangeFound = true;
0167     }
0168
0169     // Finds all tracks within the track range
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;  // at least one valid hit to account for
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