]>
Commit | Line | Data |
---|---|---|
da0e9ce3 | 1 | /******************************************************************************* |
2 | * Project: BaBar detector at the SLAC PEP-II B-factory | |
3 | * Package: EvtGenBase | |
0ca57c2f | 4 | * File: $Id: EvtPdf.hh,v 1.2 2009-03-16 16:40:15 robbep Exp $ |
da0e9ce3 | 5 | * Author: Alexei Dvoretskii, dvoretsk@slac.stanford.edu, 2001-2002 |
6 | * | |
7 | * Copyright (C) 2002 Caltech | |
8 | *******************************************************************************/ | |
9 | ||
10 | /* | |
11 | * All classes are templated on the point type T | |
12 | * | |
13 | * EvtPdf: | |
14 | * | |
15 | * Probability density function defined on an interval of phase-space. | |
16 | * Integral over the interval can be calculated by Monte Carlo integration. | |
17 | * Some (but not all) PDFs are analytic in the sense that they can be integrated | |
18 | * by numeric quadrature and distributions can be generated according to them. | |
19 | * | |
20 | * EvtPdfGen: | |
21 | * | |
22 | * Generator adaptor. Can be used to generate random points | |
23 | * distributed according to the PDF for analytic PDFs. | |
24 | * | |
25 | * EvtPdfPred: | |
26 | * | |
27 | * Predicate adaptor for PDFs. Can be used for generating random points distributed | |
28 | * according to the PDF for any PDF using rejection method. (See "Numerical Recipes"). | |
29 | * | |
30 | * EvtPdfUnary: | |
31 | * | |
32 | * Adapter for generic algorithms. Evaluates the PDF and returns the value | |
33 | * | |
34 | * EvtPdfDiv: | |
35 | * | |
36 | * PDF obtained by division of one PDF by another. Because the two PDFs are | |
37 | * arbitrary this PDF is not analytic. When importance sampling is used the | |
38 | * original PDF is divided by the analytic comparison function. EvtPdfDiv is | |
39 | * used to represent the modified PDF. | |
40 | */ | |
41 | ||
42 | #ifndef EVT_PDF_HH | |
43 | #define EVT_PDF_HH | |
44 | ||
45 | #include <assert.h> | |
46 | #include <stdio.h> | |
47 | #include "EvtGenBase/EvtValError.hh" | |
48 | #include "EvtGenBase/EvtPredGen.hh" | |
49 | #include "EvtGenBase/EvtStreamInputIterator.hh" | |
50 | #include "EvtGenBase/EvtPdfMax.hh" | |
51 | #include "EvtGenBase/EvtMacros.hh" | |
52 | #include "EvtGenBase/EvtRandom.hh" | |
53 | ||
54 | template <class T> class EvtPdfPred; | |
55 | template <class T> class EvtPdfGen; | |
56 | ||
57 | template <class T> class EvtPdf { | |
58 | public: | |
59 | ||
60 | EvtPdf() {} | |
61 | EvtPdf(const EvtPdf& other) : _itg(other._itg) {} | |
62 | virtual ~EvtPdf() {} | |
63 | virtual EvtPdf<T>* clone() const = 0; | |
64 | ||
65 | double evaluate(const T& p) const { | |
66 | if(p.isValid()) return pdf(p); | |
67 | else return 0.; | |
68 | } | |
69 | ||
70 | // Find PDF maximum. Points are sampled according to pc | |
71 | ||
72 | EvtPdfMax<T> findMax(const EvtPdf<T>& pc, int N); | |
73 | ||
74 | // Find generation efficiency. | |
75 | ||
76 | EvtValError findGenEff(const EvtPdf<T>& pc, int N, int nFindMax); | |
77 | ||
78 | // Analytic integration. Calls cascade down until an overridden | |
79 | // method is called. | |
80 | ||
81 | void setItg(EvtValError itg) {_itg = itg; } | |
82 | ||
83 | EvtValError getItg() const { | |
84 | if(!_itg.valueKnown()) _itg = compute_integral(); | |
85 | return _itg; | |
86 | } | |
87 | EvtValError getItg(int N) const { | |
88 | if(!_itg.valueKnown()) _itg = compute_integral(N); | |
89 | return _itg; | |
90 | } | |
91 | ||
92 | virtual EvtValError compute_integral() const | |
93 | //make sun happy - return something | |
94 | { printf("Analytic integration of PDF is not defined\n"); assert(0); return compute_integral();} | |
95 | virtual EvtValError compute_integral(int) const { return compute_integral(); } | |
96 | ||
97 | // Monte Carlo integration. | |
98 | ||
99 | EvtValError compute_mc_integral(const EvtPdf<T>& pc, int N); | |
100 | ||
101 | // Generation. Create predicate accept-reject generators. | |
102 | // nMax iterations will be used to find the maximum of the accept-reject predicate | |
103 | ||
104 | EvtPredGen<EvtPdfGen<T>,EvtPdfPred<T> > accRejGen(const EvtPdf<T>& pc, int nMax, double factor = 1.); | |
105 | ||
106 | virtual T randomPoint(); | |
107 | ||
108 | protected: | |
109 | ||
110 | virtual double pdf(const T&) const = 0; | |
111 | mutable EvtValError _itg; | |
112 | }; | |
113 | ||
114 | ||
115 | template <class T> class EvtPdfGen { | |
116 | public: | |
117 | typedef T result_type; | |
118 | ||
119 | EvtPdfGen() : _pdf(0) {} | |
120 | EvtPdfGen(const EvtPdfGen<T>& other) : | |
121 | _pdf(other._pdf ? other._pdf->clone() : 0) | |
122 | {} | |
123 | EvtPdfGen(const EvtPdf<T>& pdf) : | |
124 | _pdf(pdf.clone()) | |
125 | {} | |
126 | ~EvtPdfGen() { delete _pdf;} | |
127 | ||
128 | result_type operator()() {return _pdf->randomPoint();} | |
129 | ||
130 | private: | |
131 | ||
132 | EvtPdf<T>* _pdf; | |
133 | }; | |
134 | ||
135 | ||
136 | template <class T> class EvtPdfPred { | |
137 | public: | |
138 | typedef T argument_type; | |
139 | typedef bool result_type; | |
140 | ||
141 | EvtPdfPred() {} | |
142 | EvtPdfPred(const EvtPdf<T>& thePdf) : itsPdf(thePdf.clone()) {} | |
143 | EvtPdfPred(const EvtPdfPred& other) : COPY_PTR(itsPdf), COPY_MEM(itsPdfMax) {} | |
144 | ~EvtPdfPred() { delete itsPdf; } | |
145 | ||
146 | result_type operator()(argument_type p) | |
147 | { | |
148 | assert(itsPdf); | |
149 | assert(itsPdfMax.valueKnown()); | |
150 | ||
151 | double random = EvtRandom::Flat(0.,itsPdfMax.value()); | |
152 | return (random <= itsPdf->evaluate(p)); | |
153 | } | |
154 | ||
155 | EvtPdfMax<T> getMax() const { return itsPdfMax; } | |
156 | void setMax(const EvtPdfMax<T>& max) { itsPdfMax = max; } | |
157 | template <class InputIterator> void compute_max(InputIterator it, InputIterator end, | |
158 | double factor = 1.) | |
159 | { | |
160 | T p = *it++; | |
161 | itsPdfMax = EvtPdfMax<T>(p,itsPdf->evaluate(p)*factor); | |
162 | ||
163 | while(!(it == end)) { | |
164 | T p = *it++; | |
165 | double val = itsPdf->evaluate(p)*factor; | |
166 | if(val > itsPdfMax.value()) itsPdfMax = EvtPdfMax<T>(p,val); | |
167 | } | |
168 | } | |
169 | ||
170 | private: | |
171 | ||
172 | EvtPdf<T>* itsPdf; | |
173 | EvtPdfMax<T> itsPdfMax; | |
174 | }; | |
175 | ||
176 | ||
177 | template <class T> class EvtPdfUnary { | |
178 | public: | |
179 | typedef double result_type; | |
180 | typedef T argument_type; | |
181 | ||
182 | EvtPdfUnary() {} | |
183 | EvtPdfUnary(const EvtPdf<T>& thePdf) : itsPdf(thePdf.clone()) {} | |
184 | EvtPdfUnary(const EvtPdfUnary& other) : COPY_PTR(itsPdf) {} | |
185 | ~EvtPdfUnary() { delete itsPdf; } | |
186 | ||
187 | result_type operator()(argument_type p) | |
188 | { | |
189 | assert(itsPdf); | |
190 | double ret = itsPdf->evaluate(p); | |
191 | return ret; | |
192 | } | |
193 | ||
194 | private: | |
195 | ||
196 | EvtPdf<T>* itsPdf; | |
197 | }; | |
198 | ||
199 | ||
200 | template <class T> class EvtPdfDiv : public EvtPdf<T> { | |
201 | public: | |
202 | ||
203 | EvtPdfDiv() : itsNum(0), itsDen(0) {} | |
204 | EvtPdfDiv(const EvtPdf<T>& theNum, const EvtPdf<T>& theDen) | |
205 | : EvtPdf<T>(), itsNum(theNum.clone()), itsDen(theDen.clone()) | |
206 | {} | |
207 | EvtPdfDiv(const EvtPdfDiv<T>& other) | |
208 | : EvtPdf<T>(other), COPY_PTR(itsNum), COPY_PTR(itsDen) | |
209 | {} | |
210 | virtual ~EvtPdfDiv() { delete itsNum; delete itsDen; } | |
211 | virtual EvtPdf<T>* clone() const | |
212 | { return new EvtPdfDiv(*this); } | |
213 | ||
214 | virtual double pdf(const T& p) const | |
215 | { | |
216 | double num = itsNum->evaluate(p); | |
217 | double den = itsDen->evaluate(p); | |
218 | assert(den != 0); | |
219 | return num/den; | |
220 | } | |
221 | ||
222 | private: | |
223 | ||
224 | EvtPdf<T>* itsNum; // numerator | |
225 | EvtPdf<T>* itsDen; // denominator | |
226 | }; | |
227 | ||
228 | ||
229 | template <class T> | |
230 | EvtPdfMax<T> EvtPdf<T>::findMax(const EvtPdf<T>& pc, int N) | |
231 | { | |
232 | EvtPdfPred<T> pred(*this); | |
233 | EvtPdfGen<T> gen(pc); | |
234 | pred.compute_max(iter(gen,N),iter(gen)); | |
235 | EvtPdfMax<T> p = pred.getMax(); | |
236 | return p; | |
237 | } | |
238 | ||
239 | ||
240 | template <class T> | |
241 | EvtValError EvtPdf<T>::findGenEff(const EvtPdf<T>& pc, int N, int nFindMax) | |
242 | { | |
243 | assert(N > 0 || nFindMax > 0); | |
244 | EvtPredGen<EvtPdfGen<T>,EvtPdfPred<T> > gen = accRejGen(pc,nFindMax); | |
245 | int i; | |
246 | for(i=0;i<N;i++) gen(); | |
247 | double eff = double(gen.getPassed())/double(gen.getTried()); | |
248 | double err = sqrt(double(gen.getPassed()))/double(gen.getTried()); | |
249 | return EvtValError(eff,err); | |
250 | } | |
251 | ||
252 | template <class T> | |
253 | EvtValError EvtPdf<T>::compute_mc_integral(const EvtPdf<T>& pc, int N) | |
254 | { | |
255 | assert(N > 0); | |
256 | ||
257 | EvtValError otherItg = pc.getItg(); | |
258 | EvtPdfDiv<T> pdfdiv(*this,pc); | |
259 | EvtPdfUnary<T> unary(pdfdiv); | |
260 | ||
261 | EvtPdfGen<T> gen(pc); | |
262 | EvtStreamInputIterator<T> begin = iter(gen,N); | |
263 | EvtStreamInputIterator<T> end; | |
264 | ||
265 | double sum = 0.; | |
266 | double sum2 = 0.; | |
267 | while(!(begin == end)) { | |
268 | ||
269 | double value = pdfdiv.evaluate(*begin++); | |
270 | sum += value; | |
271 | sum2 += value*value; | |
272 | } | |
273 | ||
274 | EvtValError x; | |
275 | if(N > 0) { | |
276 | double av = sum/((double) N); | |
277 | if(N > 1) { | |
278 | double dev2 = (sum2 - av*av*N)/((double) (N - 1)); | |
279 | // Due to numerical precision dev2 may sometimes be negative | |
280 | if(dev2 < 0.) dev2 = 0.; | |
281 | double error = sqrt(dev2/((double) N)); | |
282 | x = EvtValError(av,error); | |
283 | } | |
284 | else x = EvtValError(av); | |
285 | } | |
286 | _itg = x * pc.getItg(); | |
287 | return _itg; | |
288 | } | |
289 | ||
290 | template <class T> | |
291 | T EvtPdf<T>::randomPoint() | |
292 | { | |
293 | printf("Function defined for analytic PDFs only\n"); | |
294 | assert(0); | |
295 | T temp; | |
296 | return temp; | |
297 | } | |
298 | ||
299 | template <class T> | |
300 | EvtPredGen<EvtPdfGen<T>,EvtPdfPred<T> > | |
301 | EvtPdf<T>::accRejGen(const EvtPdf<T>& pc, int nMax, double factor) | |
302 | { | |
303 | EvtPdfGen<T> gen(pc); | |
304 | EvtPdfDiv<T> pdfdiv(*this,pc); | |
305 | EvtPdfPred<T> pred(pdfdiv); | |
306 | pred.compute_max(iter(gen,nMax),iter(gen),factor); | |
307 | return EvtPredGen<EvtPdfGen<T>,EvtPdfPred<T> >(gen,pred); | |
308 | } | |
309 | ||
310 | #endif | |
311 | ||
312 | ||
313 | ||
314 | ||
315 |