SiStripRecHitConverter/src/CrosstalkInversion.cc

0001 #include "RecoLocalTracker/SiStripRecHitConverter/interface/CrosstalkInversion.h"
0002 #include "DataFormats/SiStripCluster/interface/SiStripCluster.h"
0003 #include "FWCore/Utilities/interface/isFinite.h"
0004
0005 namespace reco {
0006
0007   std::vector<stats_t<float> > InverseCrosstalkMatrix::unfold(const SiStripCluster& clus, const float x) {
0008     auto const& q = clus.amplitudes();
0009     const stats_t<float> suppressed(-5, 100);
0010     const stats_t<float> saturated(254, 400);
0011 #define STATS(value) ((value < 254) ? stats_t<float>(value) : saturated)
0012
0013     const unsigned N = q.size();
0014     std::vector<stats_t<float> > Q(N + 2, stats_t<float>(0));
0015     Q[0] = Q[N + 1] = suppressed;
0016
0017     if (N == 1)  //optimize N==1
0018       Q[1] = STATS(q[0]) / (1 - 2 * x);
0019     else if (N == 2) {  //optimize N==2
0020       const double A = 1 - 2 * x;
0021       const stats_t<float> q0 = STATS(q[0]);
0022       const stats_t<float> q1 = STATS(q[1]);
0023       Q[1] = (A * q0 - x * q1) / (A * A - x * x);
0024       Q[2] = (A * q1 - x * q0) / (A * A - x * x);
0025     } else {  //general case
0026       const InverseCrosstalkMatrix inverse(N, x);
0027       for (unsigned i = 0; i < (N + 1) / 2; i++) {
0028         for (unsigned j = i; j < N - i; j++) {
0029           const float Cij = inverse(i + 1, j + 1);
0030           Q[i + 1] += Cij * STATS(q[j]);
0031           if (i != j)
0032             Q[j + 1] += Cij * STATS(q[i]);
0033           if (N != i + j + 1) {
0034             Q[N - i] += Cij * STATS(q[N - j - 1]);
0035             if (i != j)
0036               Q[N - j] += Cij * STATS(q[N - i - 1]);
0037           }
0038         }
0039       }
0040     }
0041 #undef STATS
0042     return Q;
0043   }
0044
0045   InverseCrosstalkMatrix::InverseCrosstalkMatrix(const unsigned N, const float x)
0046       : N(x > 0 ? N : 0),
0047         sq(sqrt(-x * 4 + 1)),
0048         lambdaP(1 + (1 + sq) / (-x * 2)),
0049         lambdaM(1 + (1 - sq) / (-x * 2)),
0050         denominator(sq * (pow(lambdaP, N + 1) - pow(lambdaM, N + 1))) {}
0051
0052   float InverseCrosstalkMatrix::operator()(const unsigned i, const unsigned j) const {
0053     return N == 0 || edm::isNotFinite(denominator) ? i == j : i >= j ? element(i, j) : element(j, i);
0054   }
0055
0056   inline float InverseCrosstalkMatrix::element(const unsigned i, const unsigned j) const {
0057     return (pow(lambdaM, N + 1 - i) - pow(lambdaP, N + 1 - i)) * (pow(lambdaM, j) - pow(lambdaP, j)) / denominator;
0058   }
0059
0060 }  // namespace reco