TEvtGen/EvtGenBase/EvtPdf.hh

   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   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