TEvtGen/EvtGenBase/EvtPdfSum.hh

   1 /*******************************************************************************
   2  * Project: BaBar detector at the SLAC PEP-II B-factory
   3  * Package: EvtGenBase
   4  *    File: $Id: EvtPdfSum.hh,v 1.10 2009/02/19 03:22:30 ryd Exp $
   5  *  Author: Alexei Dvoretskii, dvoretsk@slac.stanford.edu, 2001-2002
   6  *
   7  * Copyright (C) 2002 Caltech
   8  *******************************************************************************/
   9
  10 // Sum of PDF functions.
  11
  12 #ifndef EVT_PDF_SUM_HH
  13 #define EVT_PDF_SUM_HH
  14
  15 #include <stdio.h>
  16 #include <vector>
  17 using std::vector;
  18 #include "EvtGenBase/EvtPdf.hh"
  19
  20 template <class T>
  21 class EvtPdfSum : public EvtPdf<T> {
  22 public:
  23
  24   EvtPdfSum() {}
  25   EvtPdfSum(const EvtPdfSum<T>& other);
  26   virtual ~EvtPdfSum();
  27   virtual EvtPdf<T>* clone() const { return new EvtPdfSum(*this); }
  28
  29
  30   // Manipulate terms and coefficients
  31
  32   void addTerm(double c,const EvtPdf<T>& pdf)
  33   { assert(c >= 0.); _c.push_back(c); _term.push_back(pdf.clone()); }
  34
  35   void addOwnedTerm(double c, EvtPdf<T>* pdf)
  36   { _c.push_back(c); _term.push_back(pdf); }
  37
  38   size_t nTerms() const { return _term.size(); }  // number of terms
  39
  40   inline double   c(int i) const { return _c[i]; }
  41   inline EvtPdf<T>* getPdf(int i) const { return _term[i]; }
  42
  43
  44   // Integrals
  45
  46   virtual EvtValError compute_integral() const;
  47   virtual EvtValError compute_integral(int N) const;
  48   virtual T randomPoint();
  49
  50 protected:
  51
  52   virtual double pdf(const T& p) const;
  53
  54   vector<double> _c;                     // coefficients
  55   vector<EvtPdf<T>*> _term;       // pointers to pdfs
  56 };
  57
  58
  59 template <class T>
  60 EvtPdfSum<T>::EvtPdfSum(const EvtPdfSum<T>& other)
  61   : EvtPdf<T>(other)
  62 {
  63   for(size_t i = 0; i < other.nTerms(); i++) {
  64     _c.push_back(other._c[i]);
  65     _term.push_back(other._term[i]->clone());
  66   }
  67 }
  68
  69 template <class T>
  70 EvtPdfSum<T>::~EvtPdfSum()
  71 {
  72   for(size_t i = 0; i < _c.size(); i++) {
  73     delete _term[i];
  74   }
  75 }
  76
  77
  78 template <class T>
  79 double EvtPdfSum<T>::pdf(const T& p) const
  80 {
  81   double ret = 0.;
  82   for(size_t i=0; i < _c.size(); i++) {
  83     ret += _c[i] * _term[i]->evaluate(p);
  84   }
  85   return ret;
  86 }
  87
  88 /*
  89  * Compute the sum integral by summing all term integrals.
  90  */
  91
  92 template <class T>
  93 EvtValError EvtPdfSum<T>::compute_integral() const
  94 {
  95   EvtValError itg(0.0,0.0);
  96   for(size_t i=0;i<nTerms();i++) {
  97     itg += _c[i]*_term[i]->getItg();
  98   }
  99   return itg;
 100 }
 101
 102 template <class T>
 103 EvtValError EvtPdfSum<T>::compute_integral(int N) const
 104 {
 105   EvtValError itg(0.0,0.0);
 106   for(size_t i=0;i<nTerms();i++) itg += _c[i]*_term[i]->getItg(N);
 107   return itg;
 108 }
 109
 110
 111 /*
 112  * Sample points randomly according to the sum of PDFs. First throw a random number uniformly
 113  * between zero and the value of the sum integral. Using this random number select one
 114  * of the PDFs. The generate a random point according to that PDF.
 115  */
 116
 117 template <class T>
 118 T EvtPdfSum<T>::randomPoint()
 119 {
 120   if(!this->_itg.valueKnown()) this->_itg = compute_integral();
 121
 122   double max = this->_itg.value();
 123   double rnd = EvtRandom::Flat(0,max);
 124
 125   double sum = 0.;
 126   size_t i;
 127   for(i = 0; i < nTerms(); i++) {
 128     double itg = _term[i]->getItg().value();
 129     sum += _c[i] * itg;
 130     if(sum > rnd) break;
 131   }
 132
 133   return _term[i]->randomPoint();
 134 }
 135
 136 #endif
 137
 138