HGCalRecProducers/src/ComputeClusterTime.cc

0001 #include "RecoLocalCalo/HGCalRecProducers/interface/ComputeClusterTime.h"
0002
0003 // functions to select the hits to compute the time of a given cluster
0004 // start with the only hits with timing information
0005 // average among the hits contained in the chosen time interval
0006
0007 // N.B. time is corrected wrt vtx-calorimeter distance
0008 // with straight line and light speed hypothesis
0009 // for charged tracks or heavy particles (longer track length or beta < 1)
0010 // need to correct the offset at analysis level
0011
0012 using namespace hgcalsimclustertime;
0013
0014 std::vector<size_t> decrease_sorted_indices(const std::vector<float>& v) {
0015   // initialize original index locations
0016   std::vector<size_t> idx(v.size());
0017   std::iota(idx.begin(), idx.end(), 0);
0018
0019   // sort indices based on comparing values in v (decreasing order)
0020   std::sort(idx.begin(), idx.end(), [&v](size_t i1, size_t i2) { return v[i1] < v[i2]; });
0021   return idx;
0022 };
0023
0024 ComputeClusterTime::ComputeClusterTime(float Xmin, float Xmax, float Cterm, float Aterm)
0025     : xMin_(Xmin), xMax_(Xmax), cTerm_(Cterm), aTerm_(Aterm) {
0026   if (xMin_ <= 0)
0027     xMin_ = 0.1;
0028 };
0029
0030 ComputeClusterTime::ComputeClusterTime() : xMin_(1.), xMax_(5.), cTerm_(0), aTerm_(0){};
0031
0032 void ComputeClusterTime::setParameters(float Xmin, float Xmax, float Cterm, float Aterm) {
0033   xMin_ = (Xmin > 0) ? Xmin : 0.1;
0034   xMax_ = Xmax;
0035   cTerm_ = Cterm;
0036   aTerm_ = Aterm;
0037   return;
0038 }
0039
0040 //time resolution parametrization
0041 float ComputeClusterTime::timeResolution(float x) {
0042   float funcVal = pow(aTerm_ / x, 2) + pow(cTerm_, 2);
0043   return sqrt(funcVal);
0044 }
0045
0046 float ComputeClusterTime::getTimeError(std::string type, float xVal) {
0047   if (type == "recHit") {
0048     //xVal is S/N
0049     //time is in ns units
0050     //std::cout << type << ", " << xVal << ", " << xMax_ << ", " < xMin_ << ", " << timeResolution(xMin_) << ", " << timeResolution(xVal) << std::endl;
0051
0052     if (xVal < xMin_)
0053       return timeResolution(xMin_);
0054     else if (xVal > xMax_)
0055       return cTerm_;
0056     else
0057       return timeResolution(xVal);
0058
0059     return -1;
0060   }
0061   return -1;
0062 }
0063
0064 //time-interval based on that ~210ps wide and with the highest number of hits
0065 //extension valid in high PU of taking smallest interval with (order of)68% of hits
0066 std::pair<float, float> ComputeClusterTime::fixSizeHighestDensity(
0067     std::vector<float>& time, std::vector<float> weight, unsigned int minNhits, float deltaT, float timeWidthBy) {
0068   if (time.size() < minNhits)
0069     return std::pair<float, float>(-99., -1.);
0070
0071   if (weight.empty())
0072     weight.resize(time.size(), 1.);
0073
0074   std::vector<float> t(time.size(), 0.);
0075   std::vector<float> w(time.size(), 0.);
0076   std::vector<size_t> sortedIndex = decrease_sorted_indices(time);
0077   for (std::size_t i = 0; i < sortedIndex.size(); ++i) {
0078     t[i] = time[sortedIndex[i]];
0079     w[i] = weight[sortedIndex[i]];
0080   }
0081
0082   int max_elements = 0;
0083   int start_el = 0;
0084   int end_el = 0;
0085   float timeW = 0.f;
0086   float tolerance = 0.05f;
0087
0088   for (auto start = t.begin(); start != t.end(); ++start) {
0089     const auto startRef = *start;
0090     int c = count_if(start, t.end(), [&](float el) { return el - startRef <= deltaT + tolerance; });
0091     if (c > max_elements) {
0092       max_elements = c;
0093       auto last_el = find_if_not(start, t.end(), [&](float el) { return el - startRef <= deltaT + tolerance; });
0094       auto valTostartDiff = *(--last_el) - startRef;
0095       if (std::abs(deltaT - valTostartDiff) < tolerance) {
0096         tolerance = std::abs(deltaT - valTostartDiff);
0097       }
0098       start_el = distance(t.begin(), start);
0099       end_el = distance(t.begin(), last_el);
0100       timeW = valTostartDiff;
0101     }
0102   }
0103
0104   // further adjust time width around the chosen one based on the hits density
0105   // proved to improve the resolution: get as many hits as possible provided they are close in time
0106   float HalfTimeDiff = timeW * timeWidthBy;
0107   float sum = 0.;
0108   float num = 0;
0109   int totSize = t.size();
0110
0111   for (int ij = 0; ij <= start_el; ++ij) {
0112     if (t[ij] > (t[start_el] - HalfTimeDiff)) {
0113       for (int kl = ij; kl < totSize; ++kl) {
0114         if (t[kl] < (t[end_el] + HalfTimeDiff)) {
0115           sum += t[kl] * w[kl];
0116           num += w[kl];
0117         } else
0118           break;
0119       }
0120       break;
0121     }
0122   }
0123
0124   if (num == 0) {
0125     return std::pair<float, float>(-99., -1.);
0126   }
0127   return std::pair<float, float>(sum / num, 1. / sqrt(num));
0128 }