Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:26:52

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   int numberOfTracks = 0;
0128   const float invResolution = 1. / resolution_;
0129   const float profileRangeMargin = sigma_ * 3.;
0130   const float profileRangeBegin = param.rangeBegin - profileRangeMargin;
0131   const float profileRangeEnd = param.rangeEnd + profileRangeMargin;
0132 
0133   std::vector<float> hitProfile((profileRangeEnd - profileRangeBegin) * invResolution + 1, 0.);
0134   // extra component to make sure that the profile drops below the threshold at range's end
0135   *hitProfile.rbegin() = -1.f;
0136 
0137   // Creates hit profile
0138   for (auto const& hit : hits) {
0139     const float center = getHitCenter(hit), rangeWidth = getHitRangeWidth(hit);
0140     std::vector<double> params{center, rangeWidth, sigma_};
0141     for (unsigned int i = 0; i < hitProfile.size(); ++i)
0142       hitProfile[i] +=
0143           pixelEfficiencyFunction_.evaluate(std::vector<double>{profileRangeBegin + i * resolution_}, params);
0144   }
0145 
0146   bool underThreshold = true;
0147   float rangeMaximum = -1.0f;
0148   bool trackRangeFound = false;
0149   int trackRangeBegin = 0, trackRangeEnd = 0;
0150 
0151   // Searches for tracks in the hit profile
0152   for (unsigned int i = 0; i < hitProfile.size(); i++) {
0153     if (hitProfile[i] > rangeMaximum)
0154       rangeMaximum = hitProfile[i];
0155 
0156     // Going above the threshold
0157     if (underThreshold && hitProfile[i] > threshold_) {
0158       underThreshold = false;
0159       trackRangeBegin = i;
0160       rangeMaximum = hitProfile[i];
0161     }
0162 
0163     // Going under the threshold
0164     else if (!underThreshold && hitProfile[i] <= threshold_) {
0165       underThreshold = true;
0166       trackRangeEnd = i;
0167       trackRangeFound = true;
0168     }
0169 
0170     // Finds all tracks within the track range
0171     if (trackRangeFound) {
0172       float trackThreshold = rangeMaximum - thresholdFromMaximum_;
0173       int trackBegin;
0174       bool underTrackThreshold = true;
0175 
0176       for (int j = trackRangeBegin; j <= trackRangeEnd; j++) {
0177         if (underTrackThreshold && hitProfile[j] > trackThreshold) {
0178           underTrackThreshold = false;
0179           trackBegin = j;
0180         } else if (!underTrackThreshold && hitProfile[j] <= trackThreshold) {
0181           underTrackThreshold = true;
0182           TRACK_TYPE track;
0183           float leftMargin = profileRangeBegin + resolution_ * trackBegin;
0184           float rightMargin = profileRangeBegin + resolution_ * j;
0185           setTrackCenter(track, 0.5f * (leftMargin + rightMargin));
0186           setTrackSigma(track, 0.5f * (rightMargin - leftMargin));
0187           result.push_back(track);
0188           numberOfTracks++;
0189         }
0190       }
0191       trackRangeFound = false;
0192     }
0193   }
0194 }
0195 
0196 template <class TRACK_TYPE, class HIT_TYPE>
0197 inline typename CTPPSTimingTrackRecognition<TRACK_TYPE, HIT_TYPE>::SpatialRange
0198 CTPPSTimingTrackRecognition<TRACK_TYPE, HIT_TYPE>::getHitSpatialRange(const HitVector& hits) {
0199   bool initialized = false;
0200   SpatialRange result = {};
0201 
0202   for (const auto& hit : hits) {
0203     const float xBegin = hit.x() - 0.5f * hit.xWidth(), xEnd = hit.x() + 0.5f * hit.xWidth();
0204     const float yBegin = hit.y() - 0.5f * hit.yWidth(), yEnd = hit.y() + 0.5f * hit.yWidth();
0205     const float zBegin = hit.z() - 0.5f * hit.zWidth(), zEnd = hit.z() + 0.5f * hit.zWidth();
0206 
0207     if (!initialized) {
0208       result.xBegin = xBegin;
0209       result.xEnd = xEnd;
0210       result.yBegin = yBegin;
0211       result.yEnd = yEnd;
0212       result.zBegin = zBegin;
0213       result.zEnd = zEnd;
0214       initialized = true;
0215       continue;
0216     }
0217     result.xBegin = std::min(result.xBegin, xBegin);
0218     result.xEnd = std::max(result.xEnd, xEnd);
0219     result.yBegin = std::min(result.yBegin, yBegin);
0220     result.yEnd = std::max(result.yEnd, yEnd);
0221     result.zBegin = std::min(result.zBegin, zBegin);
0222     result.zEnd = std::max(result.zEnd, zEnd);
0223   }
0224 
0225   return result;
0226 }
0227 
0228 template <class TRACK_TYPE, class HIT_TYPE>
0229 inline bool CTPPSTimingTrackRecognition<TRACK_TYPE, HIT_TYPE>::timeEval(const HitVector& hits,
0230                                                                         float& mean_time,
0231                                                                         float& time_sigma) const {
0232   float mean_num = 0.f, mean_denom = 0.f;
0233   bool valid_hits = false;
0234   for (const auto& hit : hits) {
0235     if (hit.tPrecision() <= 0.)
0236       continue;
0237     const float weight = std::pow(hit.tPrecision(), -2);
0238     mean_num += weight * hit.time();
0239     mean_denom += weight;
0240     valid_hits = true;  // at least one valid hit to account for
0241   }
0242   mean_time = valid_hits ? (mean_num / mean_denom) : 0.f;
0243   time_sigma = valid_hits ? std::sqrt(1.f / mean_denom) : 0.f;
0244   return valid_hits;
0245 }
0246 
0247 #endif