DeDx/interface/LikelihoodFitDeDxEstimator.h

0001 #ifndef RecoTracker_DeDx_LikelihoodFitDeDxEstimator_h
0002 #define RecoTracker_DeDx_LikelihoodFitDeDxEstimator_h
0003
0004 #include "RecoTracker/DeDx/interface/BaseDeDxEstimator.h"
0005 #include "DataFormats/TrackReco/interface/DeDxHit.h"
0006
0007 class LikelihoodFitDeDxEstimator : public BaseDeDxEstimator {
0008 public:
0009   LikelihoodFitDeDxEstimator(const edm::ParameterSet& iConfig) {}
0010
0011   std::pair<float, float> dedx(const reco::DeDxHitCollection& Hits) override {
0012     if (Hits.empty())
0013       return {0., 0.};
0014     // compute original
0015     std::array<double, 2> value;
0016     const auto& chi2 = estimate(Hits, value);
0017     // try to remove lowest dE/dx measurement
0018     const auto& n = Hits.size();
0019     if (n >= 3 && (chi2 > 1.3 * n + 4 * std::sqrt(1.3 * n))) {
0020       auto hs = Hits;
0021       hs.erase(std::min_element(hs.begin(), hs.end()));
0022       // if got better, accept
0023       std::array<double, 2> v;
0024       if (estimate(hs, v) < chi2 - 12)
0025         value = v;
0026     }
0027     return {value[0], std::sqrt(value[1])};
0028   }
0029
0030 private:
0031   void calculate_wrt_epsilon(const reco::DeDxHit&, const double&, std::array<double, 3>&);
0032   void functionEpsilon(const reco::DeDxHitCollection&, const double&, std::array<double, 3>&);
0033   double minimizeAllSaturated(const reco::DeDxHitCollection&, std::array<double, 2>&);
0034   double newtonMethodEpsilon(const reco::DeDxHitCollection&, std::array<double, 2>&);
0035   double estimate(const reco::DeDxHitCollection&, std::array<double, 2>&);
0036 };
0037
0038 /*****************************************************************************/
0039 void LikelihoodFitDeDxEstimator::calculate_wrt_epsilon(const reco::DeDxHit& h,
0040                                                        const double& epsilon,
0041                                                        std::array<double, 3>& value) {
0042   const auto& ls = h.pathLength();
0043   const auto& sn = h.error();      // energy sigma
0044   const auto y = h.charge() * ls;  // = g * y
0045   const auto sD = 2.E-3 + 0.095 * y;
0046   const auto ss = sD * sD + sn * sn;
0047   const auto s = std::sqrt(ss);
0048   const auto delta = epsilon * ls;
0049   const auto dy = delta - y;
0050   constexpr double nu(0.65);
0051
0052   // calculate derivatives with respect to delta
0053   std::array<double, 3> val{{0.}};
0054   if (h.rawDetId() == 0) {  // normal
0055     if (dy < -nu * s) {
0056       val[0] = -2. * nu * dy / s - nu * nu;
0057       val[1] = -2. * nu / s;
0058       val[2] = 0.;
0059     } else {
0060       val[0] = dy * dy / ss;
0061       val[1] = 2. * dy / ss;
0062       val[2] = 2. / ss;
0063     }
0064   } else {  // saturated
0065     if (dy < s) {
0066       val[0] = -dy / s + 1.;
0067       val[1] = -1. / s;
0068       val[2] = 0.;
0069     } else {
0070       val[0] = 0.;
0071       val[1] = 0.;
0072       val[2] = 0.;
0073     }
0074   }
0075
0076   // d/d delta -> d/d epsilon
0077   val[1] *= ls;
0078   val[2] *= ls * ls;
0079
0080   // sum
0081   for (size_t k = 0; k < value.size(); k++)
0082     value[k] += val[k];
0083 }
0084
0085 /*****************************************************************************/
0086 void LikelihoodFitDeDxEstimator::functionEpsilon(const reco::DeDxHitCollection& Hits,
0087                                                  const double& epsilon,
0088                                                  std::array<double, 3>& val) {
0089   val = {{0, 0, 0}};
0090   for (const auto& h : Hits)
0091     calculate_wrt_epsilon(h, epsilon, val);
0092 }
0093
0094 /*****************************************************************************/
0095 double LikelihoodFitDeDxEstimator::minimizeAllSaturated(const reco::DeDxHitCollection& Hits,
0096                                                         std::array<double, 2>& value) {
0097   int nStep(0);
0098   double par(3.0);  // input MeV/cm
0099
0100   std::array<double, 3> val{{0}};
0101   do {
0102     functionEpsilon(Hits, par, val);
0103     if (val[1] != 0)
0104       par += -val[0] / val[1];
0105     nStep++;
0106   } while (val[0] > 1e-3 && val[1] != 0 && nStep < 10);
0107
0108   value[0] = par * 1.1;
0109   value[1] = par * par * 0.01;
0110
0111   return val[0];
0112 }
0113
0114 /*****************************************************************************/
0115 double LikelihoodFitDeDxEstimator::newtonMethodEpsilon(const reco::DeDxHitCollection& Hits,
0116                                                        std::array<double, 2>& value) {
0117   int nStep(0);
0118   double par(3.0);  // input MeV/cm
0119   double dpar(0);
0120
0121   std::array<double, 3> val{{0}};
0122   do {
0123     functionEpsilon(Hits, par, val);
0124     if (val[2] != 0.)
0125       dpar = -val[1] / std::abs(val[2]);
0126     else
0127       dpar = 1.;  // step up, for epsilon
0128     if (par + dpar > 0)
0129       par += dpar;  // ok
0130     else
0131       par /= 2.;  // half
0132     nStep++;
0133   } while (std::abs(dpar) > 1e-3 && nStep < 50);
0134
0135   value[0] = par;
0136   value[1] = 2. / val[2];
0137
0138   return val[0];
0139 }
0140
0141 /*****************************************************************************/
0142 double LikelihoodFitDeDxEstimator::estimate(const reco::DeDxHitCollection& Hits, std::array<double, 2>& value) {
0143   // use newton method if at least one hit is not saturated
0144   for (const auto& h : Hits)
0145     if (h.rawDetId() == 0)
0146       return newtonMethodEpsilon(Hits, value);
0147   // else use minimize all saturated
0148   return minimizeAllSaturated(Hits, value);
0149 }
0150
0151 #endif