Combiner.C

Go to the documentation of this file.

#include <cmath>
#include <iterator>
#include <list>
#include <iosfwd>
// From http://www.slac.stanford.edu/~barlow/java/

struct Combiner
{
  struct Result
  {
    double fX;
    double fEl;
    double fEh;
    Result(double x=0, double el=0, double eh=0)
      : fX(x), fEl(el), fEh(eh)
    {    
    }
    virtual ~Result() {}
    double S() const
    {
      return 2 * fEl * fEh / (fEl + fEh);
    }
    double Sprime() const
    {
      return (fEh - fEl) / (fEh + fEh);
    }
    double V() const
    {
      return fEl * fEh;
    }
    double Vprime() const
    {
      return fEh - fEl;
    }
    virtual double Low() const
    {
      return fX - 3 * fEl;
    }
    virtual double High() const
    {
      return fX + 3 * fEh;
    }
  };
  
  struct Final : public Result
  {
    double fChi2;
    double fLower;
    double fUpper;
    Final(double x, double el, double eh,
    double chi2, double low, double high)
      : Result(x, el, eh),
  fChi2(chi2),
  fLower(low),
  fUpper(high)
    {}
    double Low() const
    {
      return fLower;
    }
    double High() const
    {
      return fUpper;
    }
  };    

  struct List
  {
    typedef std::list<Result> Container;
    typedef Container::iterator iterator;
    typedef Container::const_iterator const_iterator;
    
    Container fData;
    iterator begin() { return fData.begin(); }
    const_iterator begin() const { return fData.begin(); }
    iterator end() { return fData.end(); }
    const_iterator end() const { return fData.end(); }

    void Add(const Result& r)
    {
      fData.push_back(r);
    }
    void Add(double x, double el, double eh)
    {
      size_t n = fData.size();
      fData.resize(n+1);
      fData.back().fX = x;
      fData.back().fEl = el;
      fData.back().fEh = eh;
    }
  };
  typedef List::iterator iterator;
  typedef List::const_iterator const_iterator;
  
  virtual ~Combiner() {}
  
  virtual double W(const Result& r) const = 0;
  virtual double StepW(double guess, const Result& r) const = 0;
  virtual double StepOffset(double guess, const Result& r) const = 0;
  virtual double TermVar(double guess, const Result& r) const = 0;
  double Term(double guess, const Result& r) const
  {
    double var = TermVar(guess, r);
    if (var <= 0) return -1000;

    return std::pow(guess - r.fX, 2) / var;
  }
  double F(double          guess,
     double          chi2,
     const_iterator& begin,
     const_iterator& end) const
  {
    double s = -chi2;

    for (const_iterator i = begin; i != end; ++i) 
      s += Term(guess, *i);
    return s;
  }
  double FindError(unsigned short   nIter,
       const_iterator&  begin,
       const_iterator&  end,
       int              sign,
       double           best,
       double           chi2,
       double           s)
  {
    // Step size 
    double delta = 0.1 * sign * s;

    // Iterations 
    for (unsigned short i = 0; i < nIter; i++) {
      // Calculate chi^2 with current guess 
      double got = F(best + sign * s, chi2, begin, end);

      if (std::fabs(got-1) < 1e-7)
  // We're close to 1 so get out
  break;

      // The next guess' chi^2 value e
      double guess = F(best + sign * s + delta, chi2, begin, end);

      // Where to search next 
      if ((got - 1) * (guess - 1) > 0) {
  if ((got - 1) / (guess - 1) < 1)
    delta = -delta;
  else
    s += sign * delta;
  continue;
      }

      // Update error value and decrease step size 
      s     += sign * delta * (1 - got) / (guess - got);
      delta /= 2;
    }
    return s;
  }
  double FindX(unsigned short   nIter,
         const_iterator&  begin,
         const_iterator&  end,
         double           lowest,
         double           highest)
  {
    // Starting values 
    double x    = (highest+lowest)/2;
    double oldX = -1e33;

    // Do the iterations 
    for (unsigned short i = 0; i < nIter; i++) {
      double sum    = 0;
      double sumw   = 0;
      double offset = 0;

      // Loop over observations 
      for (const_iterator j = begin; j != end; ++j) {
  double w =  StepW(x, *j);
  offset   += StepOffset(x,  *j);

  sum  += j->fX * w;
  sumw += w;
      }
      x = (sum - offset) / sumw;

      if (std::fabs(x - oldX) < (highest-lowest) * 1e-6) break;
      oldX = x;
    }
    return x;
  }
  Final Calculate(const_iterator& begin,
      const_iterator& end,
      unsigned short nIter=50)
  {
    double lowest  = +1e33;
    double highest = -1e33;
    double sumLow  = 0;
    double sumHigh = 0;

    // Find boundaries and sum weights
    for (const_iterator i = begin; i != end; ++i) {
      lowest  = std::min(i->Low(),  lowest);
      highest = std::max(i->High(), highest);
      sumLow  = 1./std::pow(i->fEl, 2);
      sumHigh = 1./std::pow(i->fEh, 2);
    }
    // Summed weights 
    double sLow  = 1. / std::sqrt(sumLow);
    double sHigh = 1. / std::sqrt(sumHigh);
    
    // Now do the calculations
    double bestX    = FindX(nIter, begin, end, lowest, highest);
    double bestChi2 = F(bestX, 0, begin, end);
    double bestLow  = FindError(nIter,begin,end,-1,bestX,bestChi2,sLow);
    double bestHigh = FindError(nIter,begin,end,+1,bestX,bestChi2,sHigh);

    return Final(bestX, bestLow, bestHigh, bestChi2, lowest, highest);
  }
};

std::ostream& operator<<(std::ostream& o, const Combiner::Result& r)
{
  return o << r.fX << "\t-" << r.fEl << "\t+" << r.fEh;
}
std::ostream& operator<<(std::ostream& o, const Combiner::Final& r)
{
  return o << static_cast<const Combiner::Result&>(r)  << "\t" << r.fChi2;
}

struct LinearSigmaCombiner : public Combiner
{
  double W(const Result& r) const
  {
    double s = r.S();
    return .5 * std::pow(s + r.fX * r.Sprime(), 3) / s;
  }
  double StepW(double guess, const Result& r) const
  {
    double s = r.S();
    return s / std::pow(s+r.Sprime()*(guess - r.fX),3);
  }
  double StepOffset(double, const Result&) const
  {
    return 0;
  }
  double TermVar(double guess, const Result& r) const
  {
    return std::pow(r.S() + r.Sprime() * (guess - r.fX),2);
  }
  static double L(double guess, double x, double el, double eh)
  {
    double d  = (guess-x);
    double s  = 2 * el * eh / (el + eh);
    double sp = (eh - el) / (eh + eh);
    return std::pow(d / (s+sp*d), 2);
  }
  static double WrapL(double* xp, double* pp)
  {
    return L(xp[0], pp[0], pp[1], pp[2]);
  }
};

struct LinearVarianceCombiner : public Combiner
{
  double W(const Result& r) const
  {
    double v  = r.V();
    double vp = r.Vprime();
    return std::pow(v + r.fX * vp, 2) / (2 * v + r.fX * vp);
  }
  double StepW(double guess, const Result& r) const
  {
    double v = r.V();
    return v / std::pow(v+r.Vprime()*(guess - r.fX), 2);
  }
  double StepOffset(double guess, const Result& r) const
  {
    double vp = r.Vprime();
    return 0.5 * vp * std::pow((guess-r.fX)/(r.V()+vp*(guess-r.fX)),2);
  }
  double TermVar(double guess, const Result& r) const
  {
    return r.V() + r.Vprime() * (guess - r.fX);
  }
  static double L(double guess, double x, double el, double eh)
  {
    double d  = (guess-x);
    double v  = eh * el;
    double vp = eh - el;
    return std::pow(d,2) / (v+vp*d);
  }
  static double WrapL(double* xp, double* pp)
  {
    return L(xp[0], pp[0], pp[1], pp[2]);
  }

};

#include <TF1.h>
#include <TList.h>
#include <TLine.h>
#include <TH1.h>

struct DrawResult
{
  static TF1* MakeF(const Combiner::Result& r,
        Int_t                   j,
        Bool_t                  isSigma)
  {
    TF1* f = new TF1(Form("f%02d", j),
         isSigma ? LinearSigmaCombiner::WrapL 
         : LinearVarianceCombiner::WrapL,
         r.Low(), r.High(), 3);
    f->SetParNames("x", "#sigma^{-}", "#sigma^{+}");
    f->SetParameters(r.fX, r.fEl, r.fEh);
    f->SetLineStyle((j%3)+2);
    f->SetLineColor(kBlack);
    f->SetLineWidth(2);
    return f;
  }
  static TLine* MakeL(TF1* f)
  {
    Double_t m = f->GetParameter(0);
    TLine* l = new TLine(m-f->GetParameter(1), 1,
       m+f->GetParameter(2), 1);
    // l->SetName(Form("l%s", f->GetName()));
    l->SetLineColor(f->GetLineColor());
    l->SetLineStyle(f->GetLineStyle());
    l->SetLineWidth(f->GetLineWidth());
    return l;
  }
    
  static void Draw(Combiner::const_iterator b,
       Combiner::const_iterator e,
       Combiner::Result         r,
       Bool_t                   isSigma)
  {
    TList fs; fs.SetOwner(false);
    Int_t j = 0;
    for (Combiner::const_iterator i = b; i != e; i++) {
      TF1* f = MakeF(*i, j, isSigma);
      f->SetRange(r.Low(), r.High());
      fs.Add(f, j == 0 ? "" : "same");
      fs.Add(MakeL(f));
      j++;
    }
    TF1* fr = MakeF(r, j, isSigma);
    fr->SetLineColor(kRed+2);
    fr->SetLineStyle(1);
    fs.Add(fr);
    fs.Add(MakeL(fr));
    TIter next(&fs);
    TObject* o = 0;
    j = 0;
    while((o = next())) { o->Draw(j == 0 ? "" : "same"); j++; }
    // fs.Draw();
    static_cast<TF1*>(fs.First())->GetHistogram()
      ->GetYaxis()->SetRangeUser(-.1, 5.1);
  }
};

#include <iostream>
#include <algorithm>

void RunTest(Combiner::const_iterator b,
       Combiner::const_iterator e,
       Combiner::Final          s,
       Combiner::Final          v)
{
  std::ostream_iterator<Combiner::Result> out(std::cout, "\n");
  std::copy(b, e, out);

  LinearSigmaCombiner sc;
  Combiner::Final sf = sc.Calculate(b, e);
  Combiner::Final sd(sf.fX-s.fX, sf.fEl-s.fEl,sf.fEh-s.fEh,
         sf.fChi2-s.fChi2,0,0);
  std::cout << "Linear sigma:    " << sf  << "\n"
      << " Expected:       " << s   << "\n"
      << " Difference:     " << sd  << std::endl;
  DrawResult::Draw(b, e, sf, true);
  // return;
    
  LinearVarianceCombiner sv;
  Combiner::Final vf = sv.Calculate(b, e);
  Combiner::Final vd(vf.fX-v.fX, vf.fEl-v.fEl,vf.fEh-v.fEh,
         vf.fChi2-v.fChi2,0,0);
  std::cout << "Linear variance: " << vf  << "\n"
      << " Expected:       " << v   << "\n"
      << " Difference:     " << vd  << std::endl;

  DrawResult::Draw(b, e, vf, false);
}

       
void Test1()
{
  Combiner::List l;
  l.Add(4, 1.682, 2.346);
  l.Add(5, 1.912, 2.581);
  Combiner::Final s(4.49901, 1.33301, 1.66701, 0.115, 0, 11);
  Combiner::Final v(4.50001, 1.33701, 1.66901, 0.113, 0, 11);

  RunTest(l.begin(), l.end(), s, v);
}

void Test2()
{
  Combiner::List l;
  l.Add(6.32064,  0.567382, 0.379042);
  l.Add(6.15549,  0.159504, 0.170811);

  Combiner::Final s(6.2, 0.2, 0.2, 1, 0, 11);

  RunTest(l.begin(), l.end(), s, s);
}
void Test3()
{
  Combiner::List l;
  l.Add(6.20449,  0.451232, 0.495192);
  l.Add(6.16188,  0.165785, 0.176712);

  Combiner::Final s(6.18, 0.3, 0.3, 1, 0, 11);

  RunTest(l.begin(), l.end(), s, s);
  
}

#if 0
int
main()
{
  Test1();

  return 0;
}

#endif