Line data Source code
1 : /*******************************************************************************
2 : * Project: BaBar detector at the SLAC PEP-II B-factory
3 : * Package: EvtGenBase
4 : * File: $Id: EvtPdf.hh,v 1.2 2009-03-16 16:40:15 robbep Exp $
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 0 : EvtPdf() {}
61 0 : EvtPdf(const EvtPdf& other) : _itg(other._itg) {}
62 0 : virtual ~EvtPdf() {}
63 : virtual EvtPdf<T>* clone() const = 0;
64 :
65 : double evaluate(const T& p) const {
66 0 : if(p.isValid()) return pdf(p);
67 0 : else return 0.;
68 0 : }
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 0 : if(!_itg.valueKnown()) _itg = compute_integral();
85 0 : return _itg;
86 0 : }
87 : EvtValError getItg(int N) const {
88 0 : if(!_itg.valueKnown()) _itg = compute_integral(N);
89 0 : return _itg;
90 0 : }
91 :
92 : virtual EvtValError compute_integral() const
93 : //make sun happy - return something
94 0 : { printf("Analytic integration of PDF is not defined\n"); assert(0); return compute_integral();}
95 0 : 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 0 : _pdf(other._pdf ? other._pdf->clone() : 0)
122 0 : {}
123 : EvtPdfGen(const EvtPdf<T>& pdf) :
124 0 : _pdf(pdf.clone())
125 0 : {}
126 0 : ~EvtPdfGen() { delete _pdf;}
127 :
128 0 : 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 0 : EvtPdfPred(const EvtPdf<T>& thePdf) : itsPdf(thePdf.clone()) {}
143 0 : EvtPdfPred(const EvtPdfPred& other) : COPY_PTR(itsPdf), COPY_MEM(itsPdfMax) {}
144 0 : ~EvtPdfPred() { delete itsPdf; }
145 :
146 : result_type operator()(argument_type p)
147 : {
148 0 : assert(itsPdf);
149 0 : assert(itsPdfMax.valueKnown());
150 :
151 0 : double random = EvtRandom::Flat(0.,itsPdfMax.value());
152 0 : return (random <= itsPdf->evaluate(p));
153 : }
154 :
155 0 : EvtPdfMax<T> getMax() const { return itsPdfMax; }
156 0 : 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 0 : T p = *it++;
161 0 : itsPdfMax = EvtPdfMax<T>(p,itsPdf->evaluate(p)*factor);
162 :
163 0 : while(!(it == end)) {
164 0 : T p = *it++;
165 0 : double val = itsPdf->evaluate(p)*factor;
166 0 : if(val > itsPdfMax.value()) itsPdfMax = EvtPdfMax<T>(p,val);
167 0 : }
168 0 : }
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 0 : : EvtPdf<T>(), itsNum(theNum.clone()), itsDen(theDen.clone())
206 0 : {}
207 : EvtPdfDiv(const EvtPdfDiv<T>& other)
208 0 : : EvtPdf<T>(other), COPY_PTR(itsNum), COPY_PTR(itsDen)
209 0 : {}
210 0 : virtual ~EvtPdfDiv() { delete itsNum; delete itsDen; }
211 : virtual EvtPdf<T>* clone() const
212 0 : { return new EvtPdfDiv(*this); }
213 :
214 : virtual double pdf(const T& p) const
215 : {
216 0 : double num = itsNum->evaluate(p);
217 0 : double den = itsDen->evaluate(p);
218 0 : assert(den != 0);
219 0 : 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 0 : EvtPdfPred<T> pred(*this);
233 0 : EvtPdfGen<T> gen(pc);
234 0 : pred.compute_max(iter(gen,N),iter(gen));
235 0 : EvtPdfMax<T> p = pred.getMax();
236 : return p;
237 0 : }
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 0 : printf("Function defined for analytic PDFs only\n");
294 0 : 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 :
|
