Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:31:36

0001 #include "TrackingTools/GsfTools/interface/GSUtilities.h"
0002 
0003 #include <TMath.h>
0004 
0005 // #include <iostream>
0006 #include <cmath>
0007 
0008 #include <map>
0009 #include <functional>
0010 
0011 float GSUtilities::quantile(const float q) const {
0012   float qq = q;
0013   if (q > 1)
0014     qq = 1;
0015   else if (q < 0)
0016     qq = 0;
0017   //
0018   // mean and sigma of highest weight component
0019   //
0020   int iwMax(-1);
0021   float wMax(0.);
0022   for (unsigned i = 0; i < theNComp; i++) {
0023     if (theWeights[i] > wMax) {
0024       iwMax = i;
0025       wMax = theWeights[i];
0026     }
0027   }
0028   //
0029   // Find start values: begin with mean and
0030   // go towards f(x)=0 until two points on
0031   // either side are found (assumes monotonously
0032   // increasing function)
0033   //
0034   double x1(theParameters[iwMax]);
0035   double y1(cdf(x1) - qq);
0036   double dx = y1 > 0. ? -theErrors[iwMax] : theErrors[iwMax];
0037   double x2(x1 + dx);
0038   double y2(cdf(x2) - qq);
0039   int cnt = 0;
0040   while (y1 * y2 > 0.) {
0041     x1 = x2;
0042     y1 = y2;
0043     x2 += dx;
0044     y2 = cdf(x2) - qq;
0045     cnt++;
0046   }
0047   //   std::cout << "(" << x1 << "," << y1 << ") / ("
0048   //        << x2 << "," << y2 << ")" << std::endl;
0049   //
0050   // now use bisection to find final value
0051   //
0052   double x(0.);
0053   while (true) {
0054     // use linear interpolation
0055     x = -(x2 * y1 - x1 * y2) / (y2 - y1);
0056     double y = cdf(x) - qq;
0057     //     std::cout << x << " " << y << std::endl;
0058     if (fabs(y) < 1.e-6)
0059       break;
0060     if (y * y1 > 0.) {
0061       y1 = y;
0062       x1 = x;
0063     } else {
0064       y2 = y;
0065       x2 = x;
0066     }
0067   }
0068   return x;
0069 }
0070 
0071 float GSUtilities::errorHighestWeight() const {
0072   int iwMax(-1);
0073   float wMax(0.);
0074   for (unsigned i = 0; i < theNComp; i++) {
0075     if (theWeights[i] > wMax) {
0076       iwMax = i;
0077       wMax = theWeights[i];
0078     }
0079   }
0080   return theErrors[iwMax];
0081 }
0082 
0083 float GSUtilities::mode() const {
0084   //
0085   // start values = means of components
0086   //
0087   typedef std::multimap<double, double, std::greater<double> > StartMap;
0088   StartMap xStart;
0089   for (unsigned i = 0; i < theNComp; i++) {
0090     xStart.insert(std::pair<double, float>(pdf(theParameters[i]), theParameters[i]));
0091   }
0092   //
0093   // now try with each start value
0094   //
0095   typedef std::multimap<double, double, std::greater<double> > ResultMap;
0096   ResultMap xFound;
0097   for (StartMap::const_iterator i = xStart.begin(); i != xStart.end(); i++) {
0098     double x = findMode((*i).second);
0099     xFound.insert(std::pair<double, double>(pdf(x), x));
0100     //     std::cout << "Started at " << (*i).second
0101     //   << " , found " << x << std::endl;
0102   }
0103   //
0104   // results
0105   //
0106   //   for ( ResultMap::const_iterator i=xFound.begin();
0107   //    i!=xFound.end(); i++ ) {
0108   //     std::cout << "pdf at " << (*i).second << " = " << (*i).first << std::endl;
0109   //   }
0110   return xFound.begin()->second;
0111 }
0112 
0113 double GSUtilities::findMode(const double xStart) const {
0114   //
0115   // try with Newton
0116   //
0117   double y1(0.);
0118   double x(xStart);
0119   double y2(pdf(xStart));
0120   double yd(dpdf1(xStart));
0121   int nLoop(0);
0122   if ((y1 + y2) < 10 * DBL_MIN)
0123     return xStart;
0124   while (nLoop++ < 20 && fabs(y2 - y1) / (y2 + y1) > 1.e-6) {
0125     //     std::cout << "dy = " << y2-y1 << std::endl;
0126     double yd2 = dpdf2(x);
0127     if (fabs(yd2) < 10 * DBL_MIN)
0128       return xStart;
0129     x -= yd / dpdf2(x);
0130     yd = dpdf1(x);
0131     y1 = y2;
0132     y2 = pdf(x);
0133     //     std::cout << "New x / yd = " << x << " / " << yd << std::endl;
0134   }
0135   if (nLoop >= 20)
0136     return xStart;
0137   return x;
0138 }
0139 
0140 double GSUtilities::pdf(const double& x) const {
0141   double result(0.);
0142   for (unsigned i = 0; i < theNComp; i++)
0143     result += theWeights[i] * gauss(x, theParameters[i], theErrors[i]);
0144   return result;
0145 }
0146 
0147 double GSUtilities::cdf(const double& x) const {
0148   double result(0.);
0149   for (unsigned i = 0; i < theNComp; i++)
0150     result += theWeights[i] * gaussInt(x, theParameters[i], theErrors[i]);
0151   return result;
0152 }
0153 
0154 double GSUtilities::dpdf1(const double& x) const {
0155   double result(0.);
0156   for (unsigned i = 0; i < theNComp; i++) {
0157     double dx = (x - theParameters[i]) / theErrors[i];
0158     result += -theWeights[i] * dx / theErrors[i] * gauss(x, theParameters[i], theErrors[i]);
0159   }
0160   return result;
0161 }
0162 
0163 double GSUtilities::dpdf2(const double& x) const {
0164   double result(0.);
0165   for (unsigned i = 0; i < theNComp; i++) {
0166     double dx = (x - theParameters[i]) / theErrors[i];
0167     result += theWeights[i] / theErrors[i] / theErrors[i] * (dx * dx - 1) * gauss(x, theParameters[i], theErrors[i]);
0168   }
0169   return result;
0170 }
0171 
0172 double GSUtilities::gauss(const double& x, const double& mean, const double& sigma) const {
0173   const double fNorm(1. / sqrt(2 * TMath::Pi()));
0174   double result(0.);
0175 
0176   double d((x - mean) / sigma);
0177   if (fabs(d) < 20.)
0178     result = exp(-d * d / 2.);
0179   result *= fNorm / sigma;
0180   return result;
0181 }
0182 
0183 double GSUtilities::gaussInt(const double& x, const double& mean, const double& sigma) const {
0184   return TMath::Freq((x - mean) / sigma);
0185 }
0186 
0187 double GSUtilities::combinedMean() const {
0188   double s0(0.);
0189   double s1(0.);
0190   for (unsigned i = 0; i < theNComp; i++) {
0191     s0 += theWeights[i];
0192     s1 += theWeights[i] * theParameters[i];
0193   }
0194   return s1 / s0;
0195 }
0196 
0197 double GSUtilities::errorCombinedMean() const {
0198   double s0(0.);
0199   for (unsigned i = 0; i < theNComp; i++) {
0200     s0 += theWeights[i];
0201   }
0202   return 1. / (sqrt(s0));
0203 }
0204 
0205 float GSUtilities::maxWeight() const {
0206   // Look for the highest weight component
0207   //
0208   //int iwMax(-1);
0209   float wMax(0.);
0210   for (unsigned i = 0; i < theNComp; i++) {
0211     if (theWeights[i] > wMax) {
0212       //iwMax = i;
0213       wMax = theWeights[i];
0214     }
0215   }
0216   return wMax;
0217 }
0218 
0219 float GSUtilities::errorMode() {
0220   float mod = mode();
0221   float min = getMin(mod);
0222   float max = getMax(mod);
0223   int nBins = 1000;
0224   float dx = (max - min) / (float)nBins;
0225 
0226   float x1 = mod, x2 = mod, x1f = mod, x2f = mod;
0227   float I = 0;
0228   int cnt = 0;
0229   while (I < .68) {
0230     x1 -= dx;
0231     x2 += dx;
0232     if (pdf(x1) >= pdf(x2)) {
0233       x1f = x1;
0234       x2 -= dx;
0235     } else {
0236       x2f = x2;
0237       x1 += dx;
0238     }
0239     I = cdf(x2f) - cdf(x1f);
0240     cnt++;
0241     // for crazy pdf's return a crazy value
0242     if (cnt > 2500)
0243       return 100000.;
0244   }
0245   return (x2f - x1f) / 2;
0246 }
0247 
0248 float GSUtilities::getMin(float x) {
0249   int cnt = 0;
0250   float dx;
0251   if (fabs(x) < 2)
0252     dx = .5;
0253   else
0254     dx = fabs(x) / 10;
0255   while (cdf(x) > .1 && cnt < 1000) {
0256     x -= dx;
0257     cnt++;
0258   }
0259   return x;
0260 }
0261 
0262 float GSUtilities::getMax(float x) {
0263   int cnt = 0;
0264   float dx;
0265   if (fabs(x) < 2)
0266     dx = .5;
0267   else
0268     dx = fabs(x) / 10;
0269   while (cdf(x) < .9 && cnt < 1000) {
0270     x += dx;
0271     cnt++;
0272   }
0273   return x;
0274 }