GSUtilities.cc

Go to the documentation of this file.

#include "TrackingTools/GsfTools/interface/GSUtilities.h"

#include <TMath.h>

// #include <iostream>
#include <cmath>

#include <map>
#include <functional>

float
GSUtilities::quantile (const float q) const
{
  float qq = q;
  if (q>1) qq = 1;
  else if (q<0) qq = 0;
  //
  // mean and sigma of highest weight component
  //
  int iwMax(-1);
  float wMax(0.);
  for ( unsigned i=0; i<theNComp; i++ ) {
    if ( theWeights[i]>wMax ) {
      iwMax = i;
      wMax = theWeights[i];
    }
  }
  //
  // Find start values: begin with mean and
  // go towards f(x)=0 until two points on
  // either side are found (assumes monotonously
  // increasing function)
  //
  double x1(theParameters[iwMax]);
  double y1(cdf(x1)-qq);
  double dx = y1>0. ? -theErrors[iwMax] : theErrors[iwMax];
  double x2(x1+dx);
  double y2(cdf(x2)-qq);
  int cnt =0;
  while ( y1*y2>0. ) {
    x1 = x2;
    y1 = y2;
    x2 += dx;
    y2 = cdf(x2) - qq;
    cnt++;
  }
//   std::cout << "(" << x1 << "," << y1 << ") / (" 
//        << x2 << "," << y2 << ")" << std::endl;
  //
  // now use bisection to find final value
  //
  double x(0.);
  while ( true ) {
    // use linear interpolation
    x = -(x2*y1-x1*y2)/(y2-y1);
    double y = cdf(x) - qq;
//     std::cout << x << " " << y << std::endl;
    if ( fabs(y)<1.e-6 )  break;
    if ( y*y1>0. ) {
      y1 = y;
      x1 = x;
    }
    else {
      y2 = y;
      x2 = x;
    }
  }
  return x;
}


float
GSUtilities::errorHighestWeight() const
{
    int iwMax(-1);
    float wMax(0.);
    for ( unsigned i=0; i<theNComp; i++ ) {
        if ( theWeights[i]>wMax ) {
            iwMax = i;
            wMax = theWeights[i];
        }
    }
    return theErrors[iwMax];
}


float
GSUtilities::mode () const
{
  //
  // start values = means of components
  //
  typedef std::multimap<double, double, std::greater<double> > StartMap;
  StartMap xStart;
  for ( unsigned i=0; i<theNComp; i++ ) {
    xStart.insert(std::pair<double,float>(pdf(theParameters[i]),
                     theParameters[i]));
  }
  //
  // now try with each start value
  //
  typedef std::multimap<double, double, std::greater<double> > ResultMap;
  ResultMap xFound;
  for ( StartMap::const_iterator i=xStart.begin();
    i!=xStart.end(); i++ ) {
    double x = findMode((*i).second);
    xFound.insert(std::pair<double,double>(pdf(x),x));
//     std::cout << "Started at " << (*i).second 
//   << " , found " << x << std::endl;
  }  
  //
  // results
  //
//   for ( ResultMap::const_iterator i=xFound.begin();
//  i!=xFound.end(); i++ ) {
//     std::cout << "pdf at " << (*i).second << " = " << (*i).first << std::endl;
//   }
  return xFound.begin()->second;
}

double
GSUtilities::findMode (const double xStart) const
{
  //
  // try with Newton
  //
  double y1(0.);
  double x(xStart);
  double y2(pdf(xStart));
  double yd(dpdf1(xStart));
  int nLoop(0);
  if ( (y1+y2)<10*DBL_MIN )  return xStart;
  while ( nLoop++<20 && fabs(y2-y1)/(y2+y1)>1.e-6 ) {
//     std::cout << "dy = " << y2-y1 << std::endl;
    double yd2 = dpdf2(x);
    if ( fabs(yd2)<10*DBL_MIN )  return xStart;
    x -= yd/dpdf2(x);
    yd = dpdf1(x);
    y1 = y2;
    y2 = pdf(x);
//     std::cout << "New x / yd = " << x << " / " << yd << std::endl;
  }
  if ( nLoop >= 20 )  return xStart;
  return x;
}

double
GSUtilities::pdf (const double& x) const
{
  double result(0.);
  for ( unsigned i=0; i<theNComp; i++ )
    result += theWeights[i]*gauss(x,theParameters[i],theErrors[i]);
  return result;
}

double
GSUtilities::cdf (const double& x) const
{
  double result(0.);
  for ( unsigned i=0; i<theNComp; i++ )
    result += theWeights[i]*gaussInt(x,theParameters[i],theErrors[i]);
  return result;
}

double
GSUtilities::dpdf1 (const double& x) const
{
  double result(0.);
  for ( unsigned i=0; i<theNComp; i++ ) {
    double dx = (x-theParameters[i])/theErrors[i];
    result += -theWeights[i]*dx/theErrors[i]*
      gauss(x,theParameters[i],theErrors[i]);
  }
  return result;
}

double
GSUtilities::dpdf2 (const double& x) const
{
  double result(0.);
  for ( unsigned i=0; i<theNComp; i++ ) {
    double dx = (x-theParameters[i])/theErrors[i];
    result += theWeights[i]/theErrors[i]/theErrors[i]*
      (dx*dx-1)*gauss(x,theParameters[i],theErrors[i]);
  }
  return result;
}

double 
GSUtilities::gauss (const double& x, const double& mean,
            const double& sigma) const 
{
  const double fNorm(1./sqrt(2*TMath::Pi()));
  double result(0.);

  double d((x-mean)/sigma);
  if ( fabs(d)<20. )  result = exp(-d*d/2.);
  result *= fNorm/sigma;
  return result;
}

double 
GSUtilities::gaussInt (const double& x, const double& mean,
               const double& sigma) const 
{
  return TMath::Freq((x-mean)/sigma);
}

double
GSUtilities::combinedMean () const
{
  double s0(0.);
  double s1(0.);
  for ( unsigned i=0; i<theNComp; i++ ) {
    s0 += theWeights[i];
    s1 += theWeights[i]*theParameters[i];
  }
  return s1/s0;
}

double
GSUtilities::errorCombinedMean () const
{
  double s0(0.);
  for ( unsigned i=0; i<theNComp; i++ ) {
    s0 += theWeights[i];
  }
  return 1./(sqrt(s0));
}

float
GSUtilities::maxWeight () const
{
  // Look for the highest weight component
  //
  int iwMax(-1);
  float wMax(0.);
  for ( unsigned i=0; i<theNComp; i++ ) {
    if ( theWeights[i]>wMax ) {
      iwMax = i;
      wMax = theWeights[i];
    }
  }
  return wMax;
}


float
GSUtilities::errorMode()
{
    float mod = mode();
    float min = getMin(mod);
    float max = getMax(mod);
    int nBins = 1000;
    float dx = (max - min )/(float)nBins;

    float x1=mod, x2=mod, x1f=mod, x2f=mod;
    float I=0;
    int cnt=0;
    while (I<.68) {
        x1 -= dx;
        x2 += dx;
        if (pdf(x1)>=pdf(x2)) {
            x1f = x1;
            x2 -= dx;
        } else {
            x2f = x2;
            x1 += dx;
        }
        I = cdf(x2f) - cdf(x1f);
        cnt++;
        // for crazy pdf's return a crazy value
        if (cnt>2500) return 100000.;
    }
    return (x2f - x1f)/2;
}

float GSUtilities::getMin(float x)
{
    int cnt=0;
    float dx;
    if (fabs(x)<2) dx = .5;
    else dx = fabs(x)/10;
    while( cdf(x)>.1 && cnt<1000) {
        x -= dx;
        cnt++;
    }
    return x;
}

float GSUtilities::getMax(float x)
{
    int cnt=0;
    float dx;
    if (fabs(x)<2) dx = .5;
    else dx = fabs(x)/10;
    while( cdf(x)<.9 && cnt<1000) {
        x += dx;
        cnt++;
    }
    return x;
}