File indexing completed on 2024-04-06 12:26:29
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)
0018 Q[1] = STATS(q[0]) / (1 - 2 * x);
0019 else if (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 {
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 }