ba73bce6ce9094b465d734d0d2567f5cee56769e
[u/mrichter/AliRoot.git] / TEvtGen / EvtGenBase / EvtPdfSum.hh
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   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