Line data Source code
1 : /*******************************************************************************
2 : * Project: BaBar detector at the SLAC PEP-II B-factory
3 : * Package: EvtGenBase
4 : * File: $Id: EvtPdfSum.hh,v 1.3 2009-03-16 16:40:16 robbep 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 0 : EvtPdfSum() {}
25 0 : EvtPdfSum(const EvtPdfSum<T>& other);
26 : virtual ~EvtPdfSum();
27 0 : 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 0 : { assert(c >= 0.); _c.push_back(c); _term.push_back(pdf.clone()); }
34 :
35 : void addOwnedTerm(double c, EvtPdf<T>* pdf)
36 0 : { _c.push_back(c); _term.push_back(pdf); }
37 :
38 0 : 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 0 : : EvtPdf<T>(other)
62 0 : {
63 0 : for(size_t i = 0; i < other.nTerms(); i++) {
64 0 : _c.push_back(other._c[i]);
65 0 : _term.push_back(other._term[i]->clone());
66 : }
67 0 : }
68 :
69 : template <class T>
70 : EvtPdfSum<T>::~EvtPdfSum()
71 0 : {
72 0 : for(size_t i = 0; i < _c.size(); i++) {
73 0 : delete _term[i];
74 : }
75 0 : }
76 :
77 :
78 : template <class T>
79 : double EvtPdfSum<T>::pdf(const T& p) const
80 : {
81 : double ret = 0.;
82 0 : for(size_t i=0; i < _c.size(); i++) {
83 0 : ret += _c[i] * _term[i]->evaluate(p);
84 : }
85 0 : 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 0 : EvtValError itg(0.0,0.0);
96 0 : for(size_t i=0;i<nTerms();i++) {
97 0 : itg += _c[i]*_term[i]->getItg();
98 : }
99 : return itg;
100 0 : }
101 :
102 : template <class T>
103 : EvtValError EvtPdfSum<T>::compute_integral(int N) const
104 : {
105 0 : EvtValError itg(0.0,0.0);
106 0 : for(size_t i=0;i<nTerms();i++) itg += _c[i]*_term[i]->getItg(N);
107 : return itg;
108 0 : }
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 0 : if(!this->_itg.valueKnown()) this->_itg = compute_integral();
121 :
122 0 : double max = this->_itg.value();
123 0 : double rnd = EvtRandom::Flat(0,max);
124 :
125 : double sum = 0.;
126 : size_t i;
127 0 : for(i = 0; i < nTerms(); i++) {
128 0 : double itg = _term[i]->getItg().value();
129 0 : sum += _c[i] * itg;
130 0 : if(sum > rnd) break;
131 0 : }
132 :
133 0 : return _term[i]->randomPoint();
134 0 : }
135 :
136 : #endif
137 :
138 :
|
