BesEvtGen-00-04-08/src/EvtGen/EvtGenBase/EvtPdfSum.hh

Go to the documentation of this file.

/*******************************************************************************
 * Project: BaBar detector at the SLAC PEP-II B-factory
 * Package: EvtGenBase
 *    File: $Id: EvtPdfSum.hh,v 1.2 2007/11/20 08:36:27 pingrg Exp $
 *  Author: Alexei Dvoretskii, [email protected], 2001-2002
 *
 * Copyright (C) 2002 Caltech
 *******************************************************************************/
 
// Sum of PDF functions. 
 
#ifndef EVT_PDF_SUM_HH
#define EVT_PDF_SUM_HH
 
#include <stdio.h>
#include <vector>
using std::vector;
#include "EvtGenBase/EvtPdf.hh"
 
template <class T>
class EvtPdfSum : public EvtPdf<T> {
public:
 
  EvtPdfSum() {}
  EvtPdfSum(const EvtPdfSum<T>& other);
  virtual ~EvtPdfSum();
  virtual EvtPdf<T>* clone() const { return new EvtPdfSum(*this); }
 
 
  // Manipulate terms and coefficients
  
  void addTerm(double c,const EvtPdf<T>& pdf)
  { assert(c >= 0.); _c.push_back(c); _term.push_back(pdf.clone()); }
 
  void addOwnedTerm(double c, EvtPdf<T>* pdf)
  { _c.push_back(c); _term.push_back(pdf); }
  
  int nTerms() const { return _term.size(); }  // number of terms
  
  inline double   c(int i) const { return _c[i]; }
  inline EvtPdf<T>* getPdf(int i) const { return _term[i]; }
 
 
  // Integrals
 
  virtual EvtValError compute_integral() const;
  virtual EvtValError compute_integral(int N) const;
  virtual T randomPoint();
  
protected:
  
  virtual double pdf(const T& p) const;
  
  vector<double> _c;                     // coefficients
  vector<EvtPdf<T>*> _term;       // pointers to pdfs
  mutable EvtValError _itg;     // pingrg, 2007,Nov.19
}; 
 
 
template <class T>
EvtPdfSum<T>::EvtPdfSum(const EvtPdfSum<T>& other)
  : EvtPdf<T>(other)
{
  int i;
  for(i = 0; i < other.nTerms(); i++) {
    _c.push_back(other._c[i]);
    _term.push_back(other._term[i]->clone());
  }
}
 
template <class T>
EvtPdfSum<T>::~EvtPdfSum()
{
  int i;
  for(i = 0; i < _c.size(); i++) delete _term[i];
}
 
 
template <class T>
double EvtPdfSum<T>::pdf(const T& p) const
{
  double ret = 0.;
  unsigned i;
  for(i=0; i < _c.size(); i++) ret += _c[i] * _term[i]->evaluate(p);
  return ret;
}
 
/*
 * Compute the sum integral by summing all term integrals.
 */
 
template <class T>
EvtValError EvtPdfSum<T>::compute_integral() const 
{
  int i;
  EvtValError itg(0.0,0.0);
  for(i=0;i<nTerms();i++) itg += _c[i]*_term[i]->getItg();
  return itg;
}
 
template <class T>
EvtValError EvtPdfSum<T>::compute_integral(int N) const
{
  int i;
  EvtValError itg(0.0,0.0);
  for(i=0;i<nTerms();i++) itg += _c[i]*_term[i]->getItg(N);
  return itg;
}
 
 
/*
 * Sample points randomly according to the sum of PDFs. First throw a random number uniformly
 * between zero and the value of the sum integral. Using this random number select one
 * of the PDFs. The generate a random point according to that PDF.
 */
 
template <class T>
T EvtPdfSum<T>::randomPoint()
{
  if(!_itg.valueKnown()) _itg = compute_integral();      
  
  double max = _itg.value();
  double rnd = EvtRandom::Flat(0,max);
  
  double sum = 0.;
  int i;
  for(i = 0; i < nTerms(); i++) {
    double itg = _term[i]->getItg().value();
    sum += _c[i] * itg;
    if(sum > rnd) break;
  }
  
  return _term[i]->randomPoint(); 
}
 
#endif