GsfTools/src/GSUtilities.cc

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