AliPhysics  abafffd (abafffd)
Combiner.C
Go to the documentation of this file.
1 #include <cmath>
2 #include <iterator>
3 #include <list>
4 #include <iosfwd>
5 // From http://www.slac.stanford.edu/~barlow/java/
6 
7 struct Combiner
8 {
12  struct Result
13  {
15  double fX;
17  double fEl;
19  double fEh;
27  Result(double x=0, double el=0, double eh=0)
28  : fX(x), fEl(el), fEh(eh)
29  {
30  }
34  virtual ~Result() {}
45  double S() const
46  {
47  return 2 * fEl * fEh / (fEl + fEh);
48  }
59  double Sprime() const
60  {
61  return (fEh - fEl) / (fEh + fEh);
62  }
73  double V() const
74  {
75  return fEl * fEh;
76  }
87  double Vprime() const
88  {
89  return fEh - fEl;
90  }
96  virtual double Low() const
97  {
98  return fX - 3 * fEl;
99  }
105  virtual double High() const
106  {
107  return fX + 3 * fEh;
108  }
109  };
110 
114  struct Final : public Result
115  {
117  double fChi2;
119  double fLower;
121  double fUpper;
132  Final(double x, double el, double eh,
133  double chi2, double low, double high)
134  : Result(x, el, eh),
135  fChi2(chi2),
136  fLower(low),
137  fUpper(high)
138  {}
144  double Low() const
145  {
146  return fLower;
147  }
153  double High() const
154  {
155  return fUpper;
156  }
157  };
158 
162  struct List
163  {
165  typedef std::list<Result> Container;
167  typedef Container::iterator iterator;
169  typedef Container::const_iterator const_iterator;
170 
172  Container fData;
178  iterator begin() { return fData.begin(); }
184  const_iterator begin() const { return fData.begin(); }
190  iterator end() { return fData.end(); }
196  const_iterator end() const { return fData.end(); }
197 
203  void Add(const Result& r)
204  {
205  fData.push_back(r);
206  }
214  void Add(double x, double el, double eh)
215  {
216  size_t n = fData.size();
217  fData.resize(n+1);
218  fData.back().fX = x;
219  fData.back().fEl = el;
220  fData.back().fEh = eh;
221  }
222  };
225 
229  virtual ~Combiner() {}
230 
238  virtual double W(const Result& r) const = 0;
247  virtual double StepW(double guess, const Result& r) const = 0;
253  virtual double StepOffset(double guess, const Result& r) const = 0;
260  virtual double TermVar(double guess, const Result& r) const = 0;
276  double Term(double guess, const Result& r) const
277  {
278  double var = TermVar(guess, r);
279  if (var <= 0) return -1000;
280 
281  return std::pow(guess - r.fX, 2) / var;
282  }
294  double F(double guess,
295  double chi2,
296  const_iterator& begin,
297  const_iterator& end) const
298  {
299  double s = -chi2;
300 
301  for (const_iterator i = begin; i != end; ++i)
302  s += Term(guess, *i);
303  return s;
304  }
318  double FindError(unsigned short nIter,
319  const_iterator& begin,
320  const_iterator& end,
321  int sign,
322  double best,
323  double chi2,
324  double s)
325  {
326  // Step size
327  double delta = 0.1 * sign * s;
328 
329  // Iterations
330  for (unsigned short i = 0; i < nIter; i++) {
331  // Calculate chi^2 with current guess
332  double got = F(best + sign * s, chi2, begin, end);
333 
334  if (std::fabs(got-1) < 1e-7)
335  // We're close to 1 so get out
336  break;
337 
338  // The next guess' chi^2 value e
339  double guess = F(best + sign * s + delta, chi2, begin, end);
340 
341  // Where to search next
342  if ((got - 1) * (guess - 1) > 0) {
343  if ((got - 1) / (guess - 1) < 1)
344  delta = -delta;
345  else
346  s += sign * delta;
347  continue;
348  }
349 
350  // Update error value and decrease step size
351  s += sign * delta * (1 - got) / (guess - got);
352  delta /= 2;
353  }
354  return s;
355  }
367  double FindX(unsigned short nIter,
368  const_iterator& begin,
369  const_iterator& end,
370  double lowest,
371  double highest)
372  {
373  // Starting values
374  double x = (highest+lowest)/2;
375  double oldX = -1e33;
376 
377  // Do the iterations
378  for (unsigned short i = 0; i < nIter; i++) {
379  double sum = 0;
380  double sumw = 0;
381  double offset = 0;
382 
383  // Loop over observations
384  for (const_iterator j = begin; j != end; ++j) {
385  double w = StepW(x, *j);
386  offset += StepOffset(x, *j);
387 
388  sum += j->fX * w;
389  sumw += w;
390  }
391  x = (sum - offset) / sumw;
392 
393  if (std::fabs(x - oldX) < (highest-lowest) * 1e-6) break;
394  oldX = x;
395  }
396  return x;
397  }
407  Final Calculate(const_iterator& begin,
408  const_iterator& end,
409  unsigned short nIter=50)
410  {
411  double lowest = +1e33;
412  double highest = -1e33;
413  double sumLow = 0;
414  double sumHigh = 0;
415 
416  // Find boundaries and sum weights
417  for (const_iterator i = begin; i != end; ++i) {
418  lowest = std::min(i->Low(), lowest);
419  highest = std::max(i->High(), highest);
420  sumLow = 1./std::pow(i->fEl, 2);
421  sumHigh = 1./std::pow(i->fEh, 2);
422  }
423  // Summed weights
424  double sLow = 1. / std::sqrt(sumLow);
425  double sHigh = 1. / std::sqrt(sumHigh);
426 
427  // Now do the calculations
428  double bestX = FindX(nIter, begin, end, lowest, highest);
429  double bestChi2 = F(bestX, 0, begin, end);
430  double bestLow = FindError(nIter,begin,end,-1,bestX,bestChi2,sLow);
431  double bestHigh = FindError(nIter,begin,end,+1,bestX,bestChi2,sHigh);
432 
433  return Final(bestX, bestLow, bestHigh, bestChi2, lowest, highest);
434  }
435 };
436 
445 std::ostream& operator<<(std::ostream& o, const Combiner::Result& r)
446 {
447  return o << r.fX << "\t-" << r.fEl << "\t+" << r.fEh;
448 }
457 std::ostream& operator<<(std::ostream& o, const Combiner::Final& r)
458 {
459  return o << static_cast<const Combiner::Result&>(r) << "\t" << r.fChi2;
460 }
461 
466 {
478  double W(const Result& r) const
479  {
480  double s = r.S();
481  return .5 * std::pow(s + r.fX * r.Sprime(), 3) / s;
482  }
495  double StepW(double guess, const Result& r) const
496  {
497  double s = r.S();
498  return s / std::pow(s+r.Sprime()*(guess - r.fX),3);
499  }
505  double StepOffset(double, const Result&) const
506  {
507  return 0;
508  }
522  double TermVar(double guess, const Result& r) const
523  {
524  return std::pow(r.S() + r.Sprime() * (guess - r.fX),2);
525  }
548  static double L(double guess, double x, double el, double eh)
549  {
550  double d = (guess-x);
551  double s = 2 * el * eh / (el + eh);
552  double sp = (eh - el) / (eh + eh);
553  return std::pow(d / (s+sp*d), 2);
554  }
563  static double WrapL(double* xp, double* pp)
564  {
565  return L(xp[0], pp[0], pp[1], pp[2]);
566  }
567 };
568 
573 {
585  double W(const Result& r) const
586  {
587  double v = r.V();
588  double vp = r.Vprime();
589  return std::pow(v + r.fX * vp, 2) / (2 * v + r.fX * vp);
590  }
603  double StepW(double guess, const Result& r) const
604  {
605  double v = r.V();
606  return v / std::pow(v+r.Vprime()*(guess - r.fX), 2);
607  }
620  double StepOffset(double guess, const Result& r) const
621  {
622  double vp = r.Vprime();
623  return 0.5 * vp * std::pow((guess-r.fX)/(r.V()+vp*(guess-r.fX)),2);
624  }
638  double TermVar(double guess, const Result& r) const
639  {
640  return r.V() + r.Vprime() * (guess - r.fX);
641  }
664  static double L(double guess, double x, double el, double eh)
665  {
666  double d = (guess-x);
667  double v = eh * el;
668  double vp = eh - el;
669  return std::pow(d,2) / (v+vp*d);
670  }
679  static double WrapL(double* xp, double* pp)
680  {
681  return L(xp[0], pp[0], pp[1], pp[2]);
682  }
683 
684 };
685 
686 #include <TF1.h>
687 #include <TList.h>
688 #include <TLine.h>
689 #include <TH1.h>
690 
692 {
693  static TF1* MakeF(const Combiner::Result& r,
694  Int_t j,
695  Bool_t isSigma)
696  {
697  TF1* f = new TF1(Form("f%02d", j),
698  isSigma ? LinearSigmaCombiner::WrapL
700  r.Low(), r.High(), 3);
701  f->SetParNames("x", "#sigma^{-}", "#sigma^{+}");
702  f->SetParameters(r.fX, r.fEl, r.fEh);
703  f->SetLineStyle((j%3)+2);
704  f->SetLineColor(kBlack);
705  f->SetLineWidth(2);
706  return f;
707  }
708  static TLine* MakeL(TF1* f)
709  {
710  Double_t m = f->GetParameter(0);
711  TLine* l = new TLine(m-f->GetParameter(1), 1,
712  m+f->GetParameter(2), 1);
713  // l->SetName(Form("l%s", f->GetName()));
714  l->SetLineColor(f->GetLineColor());
715  l->SetLineStyle(f->GetLineStyle());
716  l->SetLineWidth(f->GetLineWidth());
717  return l;
718  }
719 
723  Bool_t isSigma)
724  {
725  TList fs; fs.SetOwner(false);
726  Int_t j = 0;
727  for (Combiner::const_iterator i = b; i != e; i++) {
728  TF1* f = MakeF(*i, j, isSigma);
729  f->SetRange(r.Low(), r.High());
730  fs.Add(f, j == 0 ? "" : "same");
731  fs.Add(MakeL(f));
732  j++;
733  }
734  TF1* fr = MakeF(r, j, isSigma);
735  fr->SetLineColor(kRed+2);
736  fr->SetLineStyle(1);
737  fs.Add(fr);
738  fs.Add(MakeL(fr));
739  TIter next(&fs);
740  TObject* o = 0;
741  j = 0;
742  while((o = next())) { o->Draw(j == 0 ? "" : "same"); j++; }
743  // fs.Draw();
744  static_cast<TF1*>(fs.First())->GetHistogram()
745  ->GetYaxis()->SetRangeUser(-.1, 5.1);
746  }
747 };
748 
749 #include <iostream>
750 #include <algorithm>
751 
754  Combiner::Final s,
755  Combiner::Final v)
756 {
757  std::ostream_iterator<Combiner::Result> out(std::cout, "\n");
758  std::copy(b, e, out);
759 
761  Combiner::Final sf = sc.Calculate(b, e);
762  Combiner::Final sd(sf.fX-s.fX, sf.fEl-s.fEl,sf.fEh-s.fEh,
763  sf.fChi2-s.fChi2,0,0);
764  std::cout << "Linear sigma: " << sf << "\n"
765  << " Expected: " << s << "\n"
766  << " Difference: " << sd << std::endl;
767  DrawResult::Draw(b, e, sf, true);
768  // return;
769 
771  Combiner::Final vf = sv.Calculate(b, e);
772  Combiner::Final vd(vf.fX-v.fX, vf.fEl-v.fEl,vf.fEh-v.fEh,
773  vf.fChi2-v.fChi2,0,0);
774  std::cout << "Linear variance: " << vf << "\n"
775  << " Expected: " << v << "\n"
776  << " Difference: " << vd << std::endl;
777 
778  DrawResult::Draw(b, e, vf, false);
779 }
780 
781 
782 void Test1()
783 {
784  Combiner::List l;
785  l.Add(4, 1.682, 2.346);
786  l.Add(5, 1.912, 2.581);
787  Combiner::Final s(4.49901, 1.33301, 1.66701, 0.115, 0, 11);
788  Combiner::Final v(4.50001, 1.33701, 1.66901, 0.113, 0, 11);
789 
790  RunTest(l.begin(), l.end(), s, v);
791 }
792 
793 void Test2()
794 {
795  Combiner::List l;
796  l.Add(6.32064, 0.567382, 0.379042);
797  l.Add(6.15549, 0.159504, 0.170811);
798 
799  Combiner::Final s(6.2, 0.2, 0.2, 1, 0, 11);
800 
801  RunTest(l.begin(), l.end(), s, s);
802 }
803 void Test3()
804 {
805  Combiner::List l;
806  l.Add(6.20449, 0.451232, 0.495192);
807  l.Add(6.16188, 0.165785, 0.176712);
808 
809  Combiner::Final s(6.18, 0.3, 0.3, 1, 0, 11);
810 
811  RunTest(l.begin(), l.end(), s, s);
812 
813 }
814 
815 #if 0
816 int
817 main()
818 {
819  Test1();
820 
821  return 0;
822 }
823 
824 #endif
static double WrapL(double *xp, double *pp)
Definition: Combiner.C:563
double V() const
Definition: Combiner.C:73
double Double_t
Definition: External.C:58
List::iterator iterator
Definition: Combiner.C:223
iterator begin()
Definition: Combiner.C:178
double S() const
Definition: Combiner.C:45
double Sprime() const
Definition: Combiner.C:59
double W(const Result &r) const
Definition: Combiner.C:585
void Test3()
Definition: Combiner.C:803
double Term(double guess, const Result &r) const
Definition: Combiner.C:276
Final Calculate(const_iterator &begin, const_iterator &end, unsigned short nIter=50)
Definition: Combiner.C:407
Final(double x, double el, double eh, double chi2, double low, double high)
Definition: Combiner.C:132
double TermVar(double guess, const Result &r) const
Definition: Combiner.C:522
void RunTest(Combiner::const_iterator b, Combiner::const_iterator e, Combiner::Final s, Combiner::Final v)
Definition: Combiner.C:752
Container fData
Definition: Combiner.C:172
virtual double W(const Result &r) const =0
double F(double guess, double chi2, const_iterator &begin, const_iterator &end) const
Definition: Combiner.C:294
virtual double StepOffset(double guess, const Result &r) const =0
void Test1()
Definition: Combiner.C:782
double StepOffset(double guess, const Result &r) const
Definition: Combiner.C:620
static double WrapL(double *xp, double *pp)
Definition: Combiner.C:679
List::const_iterator const_iterator
Definition: Combiner.C:224
Container::iterator iterator
Definition: Combiner.C:167
static TF1 * MakeF(const Combiner::Result &r, Int_t j, Bool_t isSigma)
Definition: Combiner.C:693
double W(const Result &r) const
Definition: Combiner.C:478
int Int_t
Definition: External.C:63
static double L(double guess, double x, double el, double eh)
Definition: Combiner.C:664
void Add(double x, double el, double eh)
Definition: Combiner.C:214
std::list< Result > Container
Definition: Combiner.C:165
void Test2()
Definition: Combiner.C:793
virtual double TermVar(double guess, const Result &r) const =0
double High() const
Definition: Combiner.C:153
std::ostream & operator<<(std::ostream &o, const Combiner::Result &r)
Definition: Combiner.C:445
static void Draw(Combiner::const_iterator b, Combiner::const_iterator e, Combiner::Result r, Bool_t isSigma)
Definition: Combiner.C:720
virtual double High() const
Definition: Combiner.C:105
iterator end()
Definition: Combiner.C:190
double fLower
Definition: Combiner.C:119
double fChi2
Definition: Combiner.C:117
double StepW(double guess, const Result &r) const
Definition: Combiner.C:603
double Low() const
Definition: Combiner.C:144
virtual ~Combiner()
Definition: Combiner.C:229
virtual double Low() const
Definition: Combiner.C:96
static double L(double guess, double x, double el, double eh)
Definition: Combiner.C:548
double StepOffset(double, const Result &) const
Definition: Combiner.C:505
const_iterator begin() const
Definition: Combiner.C:184
double FindX(unsigned short nIter, const_iterator &begin, const_iterator &end, double lowest, double highest)
Definition: Combiner.C:367
double TermVar(double guess, const Result &r) const
Definition: Combiner.C:638
void Final(Bool_t reweighed=false)
Definition: Final.C:2
double StepW(double guess, const Result &r) const
Definition: Combiner.C:495
static TLine * MakeL(TF1 *f)
Definition: Combiner.C:708
int main(int argc, char **argv)
Definition: trainMain.cxx:153
const_iterator end() const
Definition: Combiner.C:196
bool Bool_t
Definition: External.C:53
Container::const_iterator const_iterator
Definition: Combiner.C:169
double fUpper
Definition: Combiner.C:121
double Vprime() const
Definition: Combiner.C:87
double FindError(unsigned short nIter, const_iterator &begin, const_iterator &end, int sign, double best, double chi2, double s)
Definition: Combiner.C:318
virtual ~Result()
Definition: Combiner.C:34
Result(double x=0, double el=0, double eh=0)
Definition: Combiner.C:27
virtual double StepW(double guess, const Result &r) const =0
void Add(const Result &r)
Definition: Combiner.C:203