00001 /******************************************************************************* 00002 * Project: BaBar detector at the SLAC PEP-II B-factory 00003 * Package: EvtGenBase 00004 * File: $Id: EvtPdfSum.hh,v 1.2 2007/11/20 08:36:27 pingrg Exp $ 00005 * Author: Alexei Dvoretskii, dvoretsk@slac.stanford.edu, 2001-2002 00006 * 00007 * Copyright (C) 2002 Caltech 00008 *******************************************************************************/ 00009 00010 // Sum of PDF functions. 00011 00012 #ifndef EVT_PDF_SUM_HH 00013 #define EVT_PDF_SUM_HH 00014 00015 #include <stdio.h> 00016 #include <vector> 00017 using std::vector; 00018 #include "EvtGenBase/EvtPdf.hh" 00019 00020 template <class T> 00021 class EvtPdfSum : public EvtPdf<T> { 00022 public: 00023 00024 EvtPdfSum() {} 00025 EvtPdfSum(const EvtPdfSum<T>& other); 00026 virtual ~EvtPdfSum(); 00027 virtual EvtPdf<T>* clone() const { return new EvtPdfSum(*this); } 00028 00029 00030 // Manipulate terms and coefficients 00031 00032 void addTerm(double c,const EvtPdf<T>& pdf) 00033 { assert(c >= 0.); _c.push_back(c); _term.push_back(pdf.clone()); } 00034 00035 void addOwnedTerm(double c, EvtPdf<T>* pdf) 00036 { _c.push_back(c); _term.push_back(pdf); } 00037 00038 int nTerms() const { return _term.size(); } // number of terms 00039 00040 inline double c(int i) const { return _c[i]; } 00041 inline EvtPdf<T>* getPdf(int i) const { return _term[i]; } 00042 00043 00044 // Integrals 00045 00046 virtual EvtValError compute_integral() const; 00047 virtual EvtValError compute_integral(int N) const; 00048 virtual T randomPoint(); 00049 00050 protected: 00051 00052 virtual double pdf(const T& p) const; 00053 00054 vector<double> _c; // coefficients 00055 vector<EvtPdf<T>*> _term; // pointers to pdfs 00056 mutable EvtValError _itg; // pingrg, 2007,Nov.19 00057 }; 00058 00059 00060 template <class T> 00061 EvtPdfSum<T>::EvtPdfSum(const EvtPdfSum<T>& other) 00062 : EvtPdf<T>(other) 00063 { 00064 int i; 00065 for(i = 0; i < other.nTerms(); i++) { 00066 _c.push_back(other._c[i]); 00067 _term.push_back(other._term[i]->clone()); 00068 } 00069 } 00070 00071 template <class T> 00072 EvtPdfSum<T>::~EvtPdfSum() 00073 { 00074 int i; 00075 for(i = 0; i < _c.size(); i++) delete _term[i]; 00076 } 00077 00078 00079 template <class T> 00080 double EvtPdfSum<T>::pdf(const T& p) const 00081 { 00082 double ret = 0.; 00083 unsigned i; 00084 for(i=0; i < _c.size(); i++) ret += _c[i] * _term[i]->evaluate(p); 00085 return ret; 00086 } 00087 00088 /* 00089 * Compute the sum integral by summing all term integrals. 00090 */ 00091 00092 template <class T> 00093 EvtValError EvtPdfSum<T>::compute_integral() const 00094 { 00095 int i; 00096 EvtValError itg(0.0,0.0); 00097 for(i=0;i<nTerms();i++) itg += _c[i]*_term[i]->getItg(); 00098 return itg; 00099 } 00100 00101 template <class T> 00102 EvtValError EvtPdfSum<T>::compute_integral(int N) const 00103 { 00104 int i; 00105 EvtValError itg(0.0,0.0); 00106 for(i=0;i<nTerms();i++) itg += _c[i]*_term[i]->getItg(N); 00107 return itg; 00108 } 00109 00110 00111 /* 00112 * Sample points randomly according to the sum of PDFs. First throw a random number uniformly 00113 * between zero and the value of the sum integral. Using this random number select one 00114 * of the PDFs. The generate a random point according to that PDF. 00115 */ 00116 00117 template <class T> 00118 T EvtPdfSum<T>::randomPoint() 00119 { 00120 if(!_itg.valueKnown()) _itg = compute_integral(); 00121 00122 double max = _itg.value(); 00123 double rnd = EvtRandom::Flat(0,max); 00124 00125 double sum = 0.; 00126 int i; 00127 for(i = 0; i < nTerms(); i++) { 00128 double itg = _term[i]->getItg().value(); 00129 sum += _c[i] * itg; 00130 if(sum > rnd) break; 00131 } 00132 00133 return _term[i]->randomPoint(); 00134 } 00135 00136 #endif 00137 00138