]>
Commit | Line | Data |
---|---|---|
da0e9ce3 | 1 | /******************************************************************************* |
2 | * Project: BaBar detector at the SLAC PEP-II B-factory | |
3 | * Package: EvtGenBase | |
0ca57c2f | 4 | * File: $Id: EvtPdfSum.hh,v 1.3 2009-03-16 16:40:16 robbep Exp $ |
da0e9ce3 | 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 |